Box-Muller Verfahren zur Erzeugung Gaußverteilter Zufallswerte

[p. specht]

Abt. Box-Muller Verfahren zur Erzeugung Gaußverteilter Zufallswerte
------------------------------------------------------------------------------
Zufallsgeneratoren für Normalverteilung (Stichwort "Gaußsche Glockenkurve") gibts einige. Ein langsameres, aber dafür qualitativ recht gutes Verfahren haben Box und Muller 1958 vorgestellt.

WindowTitle "Box-Muller: Normalverteilte Zufallszahlen, zB. mit µ=0, "+chr$(240)+"=1"
'(D)Jan.2012 P.Specht für Paule´s PC-Forum, Demoware ohne jegliche Gewähr!
'{ Win7 Printbereich einrichten. Wer kein Win7 hat: CLS reicht!
Font 2:randomize:cls rnd(8^8)
var xr%=(%WinRight- %WinLeft-width(%Hwnd))+1
var yr%=(%WinBottom-%WinTop-height(%Hwnd))+1
var u%=640:var v%=480:var xf!=0.797:var yf!=0.72 
Window (u%-xr%)*xf!+xr%,(v%-yr%)*yf!+yr%
SetWindowPos %HWnd=(%maxx-u%)/2,((%maxy-v%)/2) - u%,v%
'}
set("decimals",16)
var pi2!= 6.283185307179586476925286766559
var fa! = 1/sqrt(pi2!)
declare r1!,r2!,z1!,z2!,s2!,s3!,n&
declare mw!,vz!,sa!,my$,my!,si5!,sigma!,s$

print " Wieviele Testzahlen werden benötigt?: ";:input n&:case n&=0:n&=10000
print " Wo liegt der reale Mittelwert?: ";:input my$:my!=val(my$)
print "\n In die sog. Sigmagrenzen fallen die folg.%-Sätze der Werte:"
print " 1"+chr$(240)+":68.27% 2:95.45% 3:99.73% 4:99.9937% 5:99.9999%"
print "\n Differenz Mittelpunkt zum Wert, der am weitesten weg liegt,"
print " [kein unplausibler Ausreisser!]  (5"+chr$(240)+") ?: ";
input si5!
sigma!=abs(si5!)/5
s$="":clearclip
r1!=rnd() ' gleichverteilt
Whileloop n&
r2!=rnd()
S2!=sqrt(-2*ln(r1!))
's3!=pi2!*r2!
z1! = S2! * cos(pi2!*r2!) '..cos(S3!) ' Standard-normalverteilt
'z2! = S2! * sin(S3!)
s$=s$+str$(z1!*sigma!+my!)
casenot &Loop=n&:s$=s$+","
mw!=mw!+z1!
vz!=vz!+sqr(z1!)
r1!=r2!
EndWhile
beep
putclip s$
print "\n Sie haben nun "+str$(n&)+" Testwerte in der Ablage!"
mw!=mw!/n&
vz!=(vz!-mw!*mw!)/(n&-1)
sa!=sqrt(vz!)
print "\n Kontrolle:  "
print " Mittelwert: ";mw!*sigma!+my!
'print " SP-Varianz: ";vz!
print " SP-Varianz: ";(sa!*sigma!)^2
print " Stand.Abw.: ";sa!*sigma!
WaitInput
End

[p. specht]

Eine Variante davon ist die von G. Marsaglia, der 1968 eine Art Wurfpfeil-Methode vorstellte. Nicht unbedingt schneller, aber simpler:

WindowTitle "Marsaglia: Normalverteilte Zufallszahlen"
' Keine Gewähr! (D) Demoware 2012 P.Specht für Paule´s PC Forum.
' Normalverteilte Zufallszahlen mit "Methode von Marsaglia"
' aka "Polarmethode": Erspart Auswertung trigonometrischer
' Funktionen, da bereits in Polarkoordinaten gerechnet wird:
CLS
set("decimals",16)
Randomize
Declare z1!,z2!
Whileloop 1000
NVZZ
Print z1!,z2!,
case %csrlin>22:cls
Endwhile
WaitInput
End

proc NVZZ
declare q!,p!,a1!,a2!
Repeat ' Hier der "Wurf":
  a1! = 2 * rnd() - 1
  a2! = 2 * rnd() - 1
  q! = a1!*a1! + a2!*a2! ' solange bis es im Kreis mit r=1 liegt
Until between(q!,0.000000000000001,1) ' Treffer!
  p! = Sqrt(-2 * ln(q!) / q!)
  z1! = a1! * p!  ' Global definierte Zufallsvariablen
  z2! = a2! * p!
endproc

P.S.: Das mit den Sigma-Grenzen (1 Sigma = Standardabweichung) kann man auch umkehren (ist vielleicht sogar plausibler): 50 % aller Messwerte haben eine Abweichung von höchstens 0,675 σ vom Mittelwert,
90 % von höchstens 1,645 σ, 95 % höchstens von 1,960 σ und 99 % höchstens von 2,575 σ vom Mittelwert.