#include "sopnamsp.h" #include "machdefs.h" #include "pexceptions.h" #include "perandom.h" #include "pemath.h" #include //////////////////////////////////////////////////////////////////////////// /*! \class SOPHYA::FunRan \ingroup NTools Classe for generating random variables from 1D function */ /********* Methode *********/ /*! Creator from a function. \verbatim - if pdf==true: f est une densite de probabilite (PDF) non necessairement normalisee if pdf==false: f est une fonction de distribution (DF). non necessairement normalisee - Le tirage aleatoire est fait sur un histogramme Histo(xMin,xMax,nBin) (voir convention dans Histo). - Chaque bin de l'histogramme contient la valeur de la PDF (ou de la DF) au centre du bin: h(i)=f(BinCenter(i)) - Les valeurs retournees sont les valeurs du centre des bins pour le tirage non interpole et toutes les valeurs entre [xmin,xmax] pour le tirage interpole - La pdf doit etre interpretee comme etant nulle pour des x<=xmin et x>=xmax - Dans le bin "i" entre [x1,x2[ et de centre x0, h(i)=pdf(x0). Pour le tirage interpole, la DF est approximee par un segment et pdf(x0) est l'exces de proba entre x1 et x2: bin 0 entre [xmin,BinHighEdge(0)[ : la pdf va de 0 a pdf(BinCenter(0)) bin 1 entre [BinLowEdge(1),BinHighEdge(1)[: la pdf va de pdf(BinCenter(0)) a pdf(BinCenter(1)) ... bin n-1 entre [BinLowEdge(n-1),xmax[: la pdf va de pdf(BinCenter(n-2)) a pdf(BinCenter(n-1)) \endverbatim */ FunRan::FunRan(FunRan::Func f, r_8 xMin, r_8 xMax, int_4 nBin, bool pdf) : Histo(xMin,xMax,nBin) { if(nBin<=1) throw RangeCheckError("FunRan::FunRan less than 2 bins requested"); for(int_4 i=0;i=0; i--) delete ranY[i]; delete[] ranY; delete ranX; } /********* Methode *********/ /*! Tirage avec retour du numeros de bin. */ void FunRan2D::BinRandom(int_4& x, int_4& y) { x = ranX->BinRandom(); y = ranY[x]->BinRandom(); } /********* Methode *********/ /*! Tirage avec retour abscisse et ordonnee du bin interpole. */ void FunRan2D::Random(r_8& x, r_8& y) { x = ranX->Random(); int_4 i = int_4(ceil(x)); y = ranY[i]->Random(); } ///////////////////////////////////////////////////////////////// /* **** Remarques sur complex< r_8 > ComplexGaussRan(double sig) **** --- variables gaussiennes x,y independantes x gaussien: pdf f(x) = 1/(sqrt(2Pi) Sx) exp(-(x-Mx)^2/(2 Sx^2)) y gaussien: pdf f(y) = 1/(sqrt(2Pi) Sy) exp(-(y-My)^2/(2 Sy^2)) x,y independants --> pdf f(x,y) = f(x) f(y) On a: = Integrate[x*f(x)] = Mx = Integrate[x^2*f(x)] = Mx^2 + Sx^2 --- On cherche la pdf g(r,t) du module et de la phase x = r cos(t) , y = r sin(t) r=sqrt(x^2+y^2 , t=atan2(y,x) (r,t) --> (x,y): le Jacobien = r g(r,t) = r f(x,y) = r f(x) f(y) = r/(2Pi Sx Sy) exp(-(x-Mx)^2/(2 Sx^2)) exp(-(y-My)^2/(2 Sy^2)) - Le cas general est complique (cf D.Pelat cours DEA "bruits et signaux" section 4.5) - Cas ou "Mx = My = 0" et "Sx = Sy = S" c'est la pdf du module et de la phase d'un nombre complexe dont les parties reelles et imaginaires sont independantes et sont distribuees selon des gaussiennes de variance S^2 g(r,t) = r/(2Pi S^2) exp(-r^2/(2 S^2)) La distribution de "r" est donc: g(r) = Integrate[g(r,t),{t,0,2Pi}] = r/S^2 exp(-r^2/(2 S^2)) La distribution de "t" est donc: g(t) = Integrate[g(r,t),{r,0,Infinity}] = 1 / 2Pi (distribution uniforme sur [0,2Pi[) Les variables aleatoires r,t sont independantes: g(r,t) = g(r) g(t) On a: = Integrate[r*g(r)] = sqrt(PI/2)*S = Integrate[r^2*g(r)] = 2*S^2 = Integrate[r^3*g(r)] = 3*sqrt(PI/2)*S^3 = Integrate[r^4*g(r)] = 8*S^4 - Attention: La variable complexe "c = x+iy = r*exp(i*t)" ainsi definie verifie: <|c|^2> = = = = 2 S^2 Si on veut generer une variable complexe gaussienne telle que = s^2 alors il faut prendre S = s/sqrt(2) comme argument */