Changeset 3109 in Sophya for trunk/SophyaLib/NTools/perandom.cc


Ignore:
Timestamp:
Nov 20, 2006, 2:14:05 PM (19 years ago)
Author:
cmv
Message:

RandomInterp() + documentation FunRan et Tirage gaussien cmv 20/11/06

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/NTools/perandom.cc

    r3097 r3109  
    1616
    1717/********* 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*/
     45FunRan::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
     56See FunRan::FunRan(FunRan::Func...) for further comments.
     57*/
     58FunRan::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.
     69If pdf=true, tab is a probability density fonction (not necessary normalised).
     70If pdf=false, tab is a distribution function (not necessarly normalized to 1).
    4571The return random values will be between 0 and nBin-1.
    4672See FunRan::FunRan(FunRan::Func...) for further comments.
    4773*/
    48 FunRan::FunRan(r_8 *tab, int_4 nBin)
     74FunRan::FunRan(r_8 *tab, int_4 nBin, bool pdf)
    4975: Histo(-0.5,nBin-0.5,nBin)
    5076{
    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.
     85If pdf=true, tab is a probability density fonction (not necessary normalised).
     86If pdf=false, tab is a distribution function (not necessarly normalized to 1).
    6587The content of tab is identified has the content of
    6688an Histogram define by Histo(xMin,xMax,nBin).
    6789See FunRan::FunRan(FunRan::Func...) for further comments.
    6890*/
    69 FunRan::FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax)
     91FunRan::FunRan(r_8 *tab, int_4 nBin, r_8 xMin, r_8 xMax, bool pdf)
    7092: Histo(xMin,xMax,nBin)
    7193{
    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
     102If pdf=true, h is a probability density fonction (not necessary normalised).
    87103If pdf=false, h is a distribution function (not necessarly normalized to 1).
    88104See FunRan::FunRan(FunRan::Func...) for further comments.
    89105*/
    90106FunRan::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 *********/
     115void 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
    100123 if((*this)(mBins-1)<=0.)
    101124   throw RangeCheckError("FunRan::FunRan(Histo) not a distribution function last bin is <=0");
    102 
    103125 for(int_4 i=0;i<mBins;i++) (*this)(i) /= (*this)(mBins-1);
    104126}
    … …  
    107129/*! Tirage avec retour du numero de bin entre 0 et mBins-1.
    108130It returns the first bin whose content is greater or equal
    109 to the random uniform number (in [0,1])
     131to the random uniform number (in [0,1[)
    110132*/
    111133int_4 FunRan::BinRandom()
    … …  
    123145 r_8 z=drand01();
    124146 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;}
    128148 return BinCenter(ibin);
    129149}
    130150
    131                    
     151/********* Methode *********/
     152/*! Tirage avec retour abscisse du bin interpole. */
     153r_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}
    132168
    133169////////////////////////////////////////////////////////////////////////////
    … …  
    219255/////////////////////////////////////////////////////////////////
    220256/*
    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
     260x gaussien: pdf f(x) = 1/(sqrt(2Pi) Sx) exp(-(x-Mx)^2/(2 Sx^2))
     261y gaussien: pdf f(y) = 1/(sqrt(2Pi) Sy) exp(-(y-My)^2/(2 Sy^2))
    224262x,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
     263On 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
    227272  g(r,t) = r f(x,y) = r f(x) f(y)
    228273         = r/(2Pi Sx Sy) exp(-(x-Mx)^2/(2 Sx^2)) exp(-(y-My)^2/(2 Sy^2))
     274
    229275- Le cas general est complique
    230276  (cf D.Pelat cours DEA "bruits et signaux" section 4.5)
     277
    231278- 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
    232282  g(r,t) = r/(2Pi S^2) exp(-r^2/(2 S^2))
    233283  La distribution de "r" est donc:
    … …  
    238288         = 1 / 2Pi  (distribution uniforme sur [0,2Pi[)
    239289  Les variables aleatoires r,t sont independantes:
    240     g(r,t) = g(r) g(t)
     290    g(r,t) = g(r) g(t)
     291On 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
    241297- Attention:
    242 La variable complexe "c=x+iy" ainsi definie verifie:
    243               <|c|^2> = <c c*> = <x^2+y^2> = <r^2> = 2 sig^2
     298La 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
    244300Si 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*/
Note: See TracChangeset for help on using the changeset viewer.