![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
Changeset 3109 in Sophya for trunk/SophyaLib/NTools
- Timestamp:
- Nov 20, 2006, 2:14:05 PM (19 years ago)
- Location:
- trunk/SophyaLib/NTools
- Files:
-
- 2 edited
-
perandom.cc (modified) (5 diffs)
-
perandom.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/NTools/perandom.cc
r3097 r3109 16 16 17 17 /********* Methode *********/ 18 /*! Createur. f is a probability density function (PDF). 19 Le tirage aleatoire est fait sur un histogramme 20 Histo(xMin,xMax,nBin) (voir convention dans Histo). 21 Chaque bin de l'histogramme contient la valeur de la PDF 22 au centre du bin. 23 Les valeurs retournees sont les valeurs du centre des bins. 24 Si binw est la largeur du bin, les valeurs retournees 25 vont de xmin+binw/2 a xmax-binw/2. 26 */ 27 FunRan::FunRan(FunRan::Func f, r_8 xMin, r_8 xMax, int_4 nBin) 28 : Histo(xMin,xMax,nBin) 29 { 30 if(nBin<0) 31 throw RangeCheckError("FunRan::FunRan less than 2 bins requested"); 32 33 // On cree la fonction de distribution a partir de la PDF 34 (*this)(0) = f(BinCenter(0)); 35 for(int_4 i=1;i<nBin;i++) (*this)(i) = (*this)(i-1) + f(BinCenter(i)); 36 37 if((*this)(nBin-1)<=0.) 38 throw RangeCheckError("FunRan::FunRan not a distribution function last bin is <=0"); 39 40 for(int_4 i=0;i<nBin;i++) (*this)(i) /= (*this)(nBin-1); 41 } 42 43 /********* Methode *********/ 44 /*! Createur. tab is a probability density function. 18 /*! Creator from a function. 19 \verbatim 20 - if pdf==true: f est une densite de probabilite (PDF) 21 non necessairement normalisee 22 if pdf==false: f est une fonction de distribution (DF). 23 non necessairement normalisee 24 - Le tirage aleatoire est fait sur un histogramme 25 Histo(xMin,xMax,nBin) (voir convention dans Histo). 26 - Chaque bin de l'histogramme contient la valeur de la PDF 27 (ou de la DF) au centre du bin: h(i)=f(BinCenter(i)) 28 - Les valeurs retournees sont les valeurs du centre des bins 29 pour le tirage non interpole et toutes les valeurs 30 entre [xmin,xmax] pour le tirage interpole 31 - La pdf doit etre interpretee comme etant nulle 32 pour des x<=xmin et x>=xmax 33 - Dans le bin "i" entre [x1,x2[ et de centre x0, h(i)=pdf(x0). 34 Pour le tirage interpole, la DF est approximee par un segment 35 et pdf(x0) est l'exces de proba entre x1 et x2: 36 bin 0 entre [xmin,BinHighEdge(0)[ : 37 la pdf va de 0 a pdf(BinCenter(0)) 38 bin 1 entre [BinLowEdge(1),BinHighEdge(1)[: 39 la pdf va de pdf(BinCenter(0)) a pdf(BinCenter(1)) 40 ... 41 bin n-1 entre [BinLowEdge(n-1),xmax[: 42 la pdf va de pdf(BinCenter(n-2)) a pdf(BinCenter(n-1)) 43 \endverbatim 44 */ 45 FunRan::FunRan(FunRan::Func f, r_8 xMin, r_8 xMax, int_4 nBin, bool pdf) 46 : Histo(xMin,xMax,nBin) 47 { 48 if(nBin<=1) 49 throw RangeCheckError("FunRan::FunRan less than 2 bins requested"); 50 for(int_4 i=0;i<nBin;i++) (*this)(i) = f(BinCenter(i)); 51 create_DF(pdf); 52 } 53 54 /********* Methode *********/ 55 /*! Creator from a ClassFunc 56 See FunRan::FunRan(FunRan::Func...) for further comments. 57 */ 58 FunRan::FunRan(ClassFunc& f, r_8 xMin, r_8 xMax, int_4 nBin, bool pdf) 59 : Histo(xMin,xMax,nBin) 60 { 61 if(nBin<=1) 62 throw RangeCheckError("FunRan::FunRan less than 2 bins requested"); 63 for(int_4 i=0;i<nBin;i++) (*this)(i) = f(BinCenter(i)); 64 create_DF(pdf); 65 } 66 67 /********* Methode *********/ 68 /*! Creator from an table. 69 If pdf=true, tab is a probability density fonction (not necessary normalised). 70 If pdf=false, tab is a distribution function (not necessarly normalized to 1). 45 71 The return random values will be between 0 and nBin-1. 46 72 See FunRan::FunRan(FunRan::Func...) for further comments. 47 73 */ 48 FunRan::FunRan(r_8 *tab, int_4 nBin )74 FunRan::FunRan(r_8 *tab, int_4 nBin, bool pdf) 49 75 : Histo(-0.5,nBin-0.5,nBin) 50 76 { 51 if(nBin<=0) 52 throw RangeCheckError("FunRan::FunRan no bin in Histo"); 53 54 (*this)(0) = tab[0]; 55 for(int_4 i=1;i<nBin;i++) (*this)(i) = (*this)(i-1) + tab[i]; 56 57 if((*this)(nBin-1)<=0.) 58 throw RangeCheckError("FunRan::FunRan not a distribution function last bin is <=0"); 59 60 for(int_4 i=0;i<nBin;i++) (*this)(i) /= (*this)(nBin-1); 61 } 62 63 /********* Methode *********/ 64 /*! Createur. tab is a probability density function 77 if(nBin<=1) 78 throw RangeCheckError("FunRan::FunRan less than 2 bins requested"); 79 for(int_4 i=0;i<nBin;i++) (*this)(i) = tab[i]; 80 create_DF(pdf); 81 } 82 83 /********* Methode *********/ 84 /*! Creator from an table. 85 If pdf=true, tab is a probability density fonction (not necessary normalised). 86 If pdf=false, tab is a distribution function (not necessarly normalized to 1). 65 87 The content of tab is identified has the content of 66 88 an Histogram define by Histo(xMin,xMax,nBin). 67 89 See FunRan::FunRan(FunRan::Func...) for further comments. 68 90 */ 69 FunRan::FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax )91 FunRan::FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax, bool pdf) 70 92 : Histo(xMin,xMax,nBin) 71 93 { 72 if(nBin<=0) 73 throw RangeCheckError("FunRan::FunRan no bin in Histo"); 74 75 (*this)(0) = tab[0]; 76 for(int_4 i=1;i<nBin;i++) (*this)(i) = (*this)(i-1) + tab[i]; 77 78 if((*this)(nBin-1)<=0.) 79 throw RangeCheckError("FunRan::FunRan not a distribution function last bin is <=0"); 80 81 for(int_4 i=0;i<nBin;i++) (*this)(i) /= (*this)(nBin-1); 82 } 83 84 /********* Methode *********/ 85 /*! Createur. 86 If pdf=true, h is a probability density fonction. 94 if(nBin<=1) 95 throw RangeCheckError("FunRan::FunRan less than 2 bins requested"); 96 for(int_4 i=0;i<nBin;i++) (*this)(i) = tab[i]; 97 create_DF(pdf); 98 } 99 100 /********* Methode *********/ 101 /*! Creator from an histogram 102 If pdf=true, h is a probability density fonction (not necessary normalised). 87 103 If pdf=false, h is a distribution function (not necessarly normalized to 1). 88 104 See FunRan::FunRan(FunRan::Func...) for further comments. 89 105 */ 90 106 FunRan::FunRan(Histo &h, bool pdf) 91 : Histo(h) 92 { 93 if(mBins<=0) 94 throw RangeCheckError("FunRan::FunRan(Histo) no bin in Histo"); 95 96 // On cree l'histo integre 97 if(pdf) 98 for(int_4 i=1;i<mBins;i++) (*this)(i) += (*this)(i-1); 99 107 : Histo(h) 108 { 109 if(mBins<=1) 110 throw RangeCheckError("FunRan::FunRan less than 2 bins requested"); 111 create_DF(pdf); 112 } 113 114 /********* Methode *********/ 115 void FunRan::create_DF(bool pdf) 116 // Creation (si necessaire) et normalisation de la DF 117 { 118 // On fabrique la FD si necessaire 119 if(pdf) 120 for(int_4 i=1;i<mBins;i++) (*this)(i) += (*this)(i-1); 121 122 // On normalise la FD 100 123 if((*this)(mBins-1)<=0.) 101 124 throw RangeCheckError("FunRan::FunRan(Histo) not a distribution function last bin is <=0"); 102 103 125 for(int_4 i=0;i<mBins;i++) (*this)(i) /= (*this)(mBins-1); 104 126 } … … 107 129 /*! Tirage avec retour du numero de bin entre 0 et mBins-1. 108 130 It returns the first bin whose content is greater or equal 109 to the random uniform number (in [0,1 ])131 to the random uniform number (in [0,1[) 110 132 */ 111 133 int_4 FunRan::BinRandom() … … 123 145 r_8 z=drand01(); 124 146 int ibin = mBins-1; 125 for(int_4 i=0;i<mBins;i++) 126 if(z<(*this)(i)) {ibin=i; break;} 127 147 for(int_4 i=0;i<mBins;i++) if(z<(*this)(i)) {ibin=i; break;} 128 148 return BinCenter(ibin); 129 149 } 130 150 131 151 /********* Methode *********/ 152 /*! Tirage avec retour abscisse du bin interpole. */ 153 r_8 FunRan::RandomInterp(void) 154 { 155 r_8 z=drand01(); 156 int ibin = mBins-1; 157 r_8 z1=0., z2; 158 for(int_4 i=0;i<mBins;i++) { 159 z2 = (*this)(i); 160 if(z<z2) {ibin=i; break;} 161 z1 = z2; 162 } 163 164 // l'algorithme garanti que "z2-z1 != 0" et "z1<z2" 165 return BinLowEdge(ibin) + (z-z1)/(z2-z1)*binWidth; 166 167 } 132 168 133 169 //////////////////////////////////////////////////////////////////////////// … … 219 255 ///////////////////////////////////////////////////////////////// 220 256 /* 221 --- Remarque sur complex< r_8 > ComplexGaussRan(double sig) 222 x = r cos(t) tire gaussien: pdf f(x) = 1/(sqrt(2Pi) Sx) exp(-(x-Mx)^2/(2 Sx^2)) 223 y = r sin(t) tire gaussien: pdf f(y) = 1/(sqrt(2Pi) Sy) exp(-(y-My)^2/(2 Sy^2)) 257 **** Remarques sur complex< r_8 > ComplexGaussRan(double sig) **** 258 259 --- variables gaussiennes x,y independantes 260 x gaussien: pdf f(x) = 1/(sqrt(2Pi) Sx) exp(-(x-Mx)^2/(2 Sx^2)) 261 y gaussien: pdf f(y) = 1/(sqrt(2Pi) Sy) exp(-(y-My)^2/(2 Sy^2)) 224 262 x,y independants --> pdf f(x,y) = f(x) f(y) 225 - On cherche la pdf g(r,t) du module et de la phase 226 (r=sqrt(x^2+y^2,t=atan2(y,x)) --> (x,y): le Jacobien = r 263 On a: 264 <x> = Integrate[x*f(x)] = Mx 265 <x^2> = Integrate[x^2*f(x)] = Mx^2 + Sx^2 266 267 --- On cherche la pdf g(r,t) du module et de la phase 268 x = r cos(t) , y = r sin(t) 269 r=sqrt(x^2+y^2 , t=atan2(y,x) 270 (r,t) --> (x,y): le Jacobien = r 271 227 272 g(r,t) = r f(x,y) = r f(x) f(y) 228 273 = r/(2Pi Sx Sy) exp(-(x-Mx)^2/(2 Sx^2)) exp(-(y-My)^2/(2 Sy^2)) 274 229 275 - Le cas general est complique 230 276 (cf D.Pelat cours DEA "bruits et signaux" section 4.5) 277 231 278 - Cas ou "Mx = My = 0" et "Sx = Sy = S" 279 c'est la pdf du module et de la phase d'un nombre complexe 280 dont les parties reelles et imaginaires sont independantes 281 et sont distribuees selon des gaussiennes de variance S^2 232 282 g(r,t) = r/(2Pi S^2) exp(-r^2/(2 S^2)) 233 283 La distribution de "r" est donc: … … 238 288 = 1 / 2Pi (distribution uniforme sur [0,2Pi[) 239 289 Les variables aleatoires r,t sont independantes: 240 g(r,t) = g(r) g(t) 290 g(r,t) = g(r) g(t) 291 On a: 292 <r> = Integrate[r*g(r)] = sqrt(PI/2)*S 293 <r^2> = Integrate[r^2*g(r)] = 2*S^2 294 <r^3> = Integrate[r^3*g(r)] = 3*sqrt(PI/2)*S^3 295 <r^4> = Integrate[r^4*g(r)] = 8*S^4 296 241 297 - Attention: 242 La variable complexe "c =x+iy" ainsi definie verifie:243 <|c|^2> = <c c*> = <x^2+y^2> = <r^2> = 2 sig^2298 La variable complexe "c = x+iy = r*exp(i*t)" ainsi definie verifie: 299 <|c|^2> = <c c*> = <x^2+y^2> = <r^2> = 2 S^2 244 300 Si on veut generer une variable complexe gaussienne telle que 245 <c c*> = s^2 alors il faut sig = s/sqrt(2) comme argument 246 */ 301 <c c*> = s^2 alors il faut prendre S = s/sqrt(2) comme argument 302 303 */ -
trunk/SophyaLib/NTools/perandom.h
r3100 r3109 12 12 #include "histos.h" 13 13 #include "srandgen.h" 14 #include "classfunc.h" 14 15 #include <complex> 15 16 … … 19 20 public: 20 21 typedef r_8 (*Func)(r_8); 21 FunRan(Func f, r_8 xMin=0.0, r_8 xMax=1.0, int_4 nBin=100); 22 FunRan(r_8 *tab, int_4 nBin); 23 FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax); 22 FunRan(ClassFunc& f, r_8 xMin=0.0, r_8 xMax=1.0, int_4 nBin=100, bool pdf=true); 23 FunRan(Func f, r_8 xMin=0.0, r_8 xMax=1.0, int_4 nBin=100, bool pdf=true); 24 FunRan(r_8 *tab, int_4 nBin, bool pdf=true); 25 FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax, bool pdf=true); 24 26 FunRan(Histo &h, bool pdf=true); 27 int_4 BinRandom(void); 25 28 r_8 Random(void); 26 int_4 BinRandom(void); 29 r_8 RandomInterp(void); 30 protected: 31 void create_DF(bool pdf); 27 32 }; 28 33
Note:
See TracChangeset
for help on using the changeset viewer.
