Metoda polarna Marsaglii

Metoda polarna Marsaglii to metoda próbkowania liczb pseudolosowych służąca do generowania pary niezależnych standardowych normalnych zmiennych losowych .

Standardowe normalne zmienne losowe są często wykorzystywane w informatyce , statystyce obliczeniowej , aw szczególności w zastosowaniach metody Monte Carlo .

Metoda biegunowa polega na wybieraniu losowych punktów ( x , y ) w kwadracie −1 < x < 1, −1 < y < 1 do

${\ Displaystyle 0 <s = x ^ {2} + y ^ {2}<1, \,}$

a następnie zwrócenie wymaganej pary normalnych zmiennych losowych jako

${\ Displaystyle x {\ sqrt {\ Frac {-2 \ ln (s)}} {s}}} \ \ \ \ y {\sqrt {\frac {-2\ln(s)}{s}}},}$

lub równoważnie,

${\ Displaystyle {\ Frac {x} {\ sqrt {s}}} {\ sqrt {-2 \ ln (s)}} \,,\ \ {\frac {y}{\sqrt {s}}}{\sqrt {-2\ln(s)}},}$

gdzie ${\ Displaystyle x/{\ sqrt {s}}}$ i ${\ Displaystyle y / {\ sqrt {s}}}$ reprezentują cosinus i sinus kąta, pod którym wektor ( x , y ) tworzy z osią x .

Podstawy teoretyczne

Podstawową teorię można podsumować w następujący sposób:

Jeśli u jest równomiernie rozłożone w przedziale 0 ≤ u < 1, to punkt (cos(2π u ), sin(2π u )) jest równomiernie rozłożony na obwodzie jednostki x ² + y ² = 1 i pomnożenie tego punktu przez niezależna zmienna losowa ρ, której rozkład to

${\ Displaystyle \ Pr (\ rho <a) = \ int _ {0} ^ {a} re ^ {- r ^ {2} / 2}\,dr}$

wygeneruje punkt

${\ Displaystyle \ lewo (\ rho \ cos (2 \ pi u), \ rho \ sin (2 \ pi u) \ prawo)}$

którego współrzędne są wspólnie rozłożone jako dwie niezależne standardowe normalne zmienne losowe.

Historia

Pomysł ten pochodzi od Laplace'a , któremu Gauss przypisuje odkrycie powyższego

${\ Displaystyle I = \ int _ {- \ infty} ^ {\ infty} e ^ {-x ^ {2}/2} \, dx}$

biorąc pierwiastek kwadratowy z

${\ Displaystyle I ^ {2} = \ int _ {- \ infty} ^ {\ infty} \ int _ {- \ infty} ^ {\ infty} e ^ {- (x ^ {2} + y ^ {2 })/2}\,dx\,dy=\int _{0}^{2\pi }\int _{0}^{\infty }re^{-r^{2}/2}\,dr \,d\theta .}$

Transformacja do współrzędnych biegunowych pokazuje, że θ jest równomiernie rozłożona (stała gęstość) od 0 do 2π, a odległość promieniowa r ma gęstość

${\ Displaystyle re ^ {- r ^ {2}/2}. \,}$

( r ² ma odpowiedni rozkład chi-kwadrat .)

Ta metoda tworzenia pary niezależnych standardowych zmiennych normalnych poprzez promieniowe rzutowanie losowego punktu na obwodzie jednostki na odległość określoną przez pierwiastek kwadratowy zmiennej chi-kwadrat-2 nazywana jest polarną metodą generowania pary normalnych zmiennych losowych ,

Względy praktyczne

Bezpośrednie zastosowanie tego pomysłu

${\ Displaystyle x = {\ sqrt {-2 \ ln (u_{1})}}\cos(2\pi u_{2}),\quad y={\sqrt {-2\ln(u_{1})}}\sin(2\pi u_{2} )}$

nazywa się transformacją Boxa-Mullera , w której zmienna chi jest zwykle generowana jako

${\ Displaystyle {\ sqrt {-2 \ ln (u_ {1})}},}$

ale ta transformacja wymaga funkcji logarytmu, pierwiastka kwadratowego, sinusa i cosinusa. W niektórych procesorach cosinus i sinus tego samego argumentu można obliczyć równolegle przy użyciu jednej instrukcji. Szczególnie w przypadku maszyn opartych na Intelu można użyć instrukcji asemblera fsincos lub instrukcji expi (dostępnej np. w D), aby obliczyć złożone

${\ Displaystyle \ nazwa operatora {wyrażenie} (z) = e ^ {iz} = \ cos (z) + i \ sin (z),\,}$

i po prostu oddziel rzeczywistą i urojoną część.

Uwaga: Aby jawnie obliczyć zespoloną postać biegunową, użyj następujących podstawień w postaci ogólnej,

Niech ${\ Displaystyle R = {\ sqrt {-2 \ ln (u_ {1})}}}$ i ${\ displaystyle z = 2 \ pi u_ {2}.}$ Następnie

${\ Displaystyle \ re ^ {iz} = {\ sqrt {-2 \ ln (u_ {1})}} e ^ {i2 \ pi u_ {2}} = {\ sqrt {-2 \ ln (u_ {1} })}}\lewo[\cos(2\pi u_{2})+i\sin(2\pi u_{2})\prawo].}$

W przeciwieństwie do tego metoda biegunowa eliminuje tutaj potrzebę obliczania cosinusa i sinusa. Zamiast tego, rozwiązując punkt na okręgu jednostkowym, te dwie funkcje można zastąpić współrzędnymi x i y znormalizowanymi do ${\ Displaystyle {\ sqrt {x ^ {2} + y ^ {2} }}}$ promień. W szczególności losowy punkt ( x , y ) wewnątrz okręgu jednostkowego jest rzutowany na obwód jednostkowy przez ustawienie ${\ displaystyle s = x ^ {2} + y ^ {2}}$ i tworząc punkt

${\ Displaystyle \ lewo ({\ Frac {x} {\ sqrt {s}}}, {\ Frac {y} {\ sqrt {s}}} \ prawej), \,}$

co jest szybszą procedurą niż obliczanie cosinusa i sinusa. Niektórzy badacze twierdzą, że instrukcja warunkowa if (do odrzucania punktu poza okręgiem jednostkowym) może spowolnić programy na nowoczesnych procesorach wyposażonych w potokowanie i przewidywanie rozgałęzień. $.$ generatora liczb losowych (tylko wygenerowanych leży wewnątrz okręgu jednostkowego)

Ten przypadkowy punkt na obwodzie jest następnie rzutowany promieniowo na wymaganą przypadkową odległość za pomocą

${\ Displaystyle {\ sqrt {-2 \ ln (s)}}, \,}$

używając tego samego s , ponieważ to s jest niezależne od przypadkowego punktu na obwodzie i samo w sobie jest równomiernie rozłożone od 0 do 1.

Realizacja

Prosta implementacja w Javie z wykorzystaniem średniej i odchylenia standardowego:

   
     

        
      
          
             
     prywatny  statyczny  podwójny  zapasowy  ;  private  static  boolean  hasSpare  =  false  ;  publiczne  statyczne  zsynchronizowane  podwójne  generowanie Gaussa  (  podwójna  średnia  ,  podwójne  stdDev  )  {  if  (  hasSpare  )  {  hasSpare  =  false  ;  zwróć  zapas  *  stdev  +  średnia  ;  }   
           
         
                  
                  
                    
               else  {  podwójne  u  ,  v  ,  s  ;  zrobić  {  u  =  matematyka  .  losowo  ()  *  2  -  1  ;  v  =  Matematyka  .  losowo  ()  *  2  -  1  ;  s  =  u  *  u  +  v  *  v  ;  }  podczas  (  s  >=  1  ||  s   0
              
            
          
               
    
 ==  );  s  =  Matematyka  .  sqrt  (  -2.0  *  Math  .  log  (  s  )  /  s  )  ;  zapas  =  v  *  s  ;  maZapas  =  true  ;  zwróć  średnią  +  stdDev  *  u  *  s  ;  }  } 

nieobsługująca wątków w C++ przy użyciu średniej i odchylenia standardowego:

     
      
        

      
          
             
      
            podwójne  generowanie Gaussa  (  podwójna  średnia  ,  podwójne  stdDev  )  {  statyczne  podwójne  zapasowe  ;  static  bool  hasSpare  =  false  ;  if  (  maZapas  )  {  maZapas  =  fałsz  ;  zwróć  zapas  *  stdev  +  średnia  ;  }  else  {  podwójne  u  ,  v  ,  s 
         
                    
                    
                    
           ;  do  {  u  =  (  rand  ()  /  ((  double  )  RAND_MAX  ))  *  2,0  -  1,0  ;  v  =  (  rand  ()  /  ((  double  )  RAND_MAX  ))  *  2,0  -  1,0  ;  s  =  u  *  u  +  v  *  v  ;  }  podczas  (  s       
              
            
          
               
    
 >=  1,0  ||  s  ==  0,0  );  s  =  sqrt  (  -2.0  *  log  (  s  )  /  s  );  zapas  =  v  *  s  ;  maZapas  =  true  ;  zwróć  średnią  +  stdDev  *  u  *  s  ;  }  } 

Implementacja std::normal_distribution w C++11 GNU GCC libstdc++ używa polarnej metody Marsaglii, cytowanej tutaj .