Sz´ am. szim. labor 2015 4. ea. T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as
Sz´am´ıt´og´epes szimul´aci´ok labor 2015 4. V´eletlen sz´ amok
HF Hivatkoz´ asok
T˝ oke Csaba BME Fizika Int´ ezet
2015. okt´ ober 7.
Sz´ am. szim. labor 2015 4. ea.
V´azlat
T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
• Egyenletes eloszl´ as´ u pszeudov´eletlen sz´amok
• V´ eletlen sz´ amok gener´ al´ asa GSL-lel
• Adott eloszl´ as szerinti v´eletlen sz´ amok gener´al´asa • V´ eletlensz´am-gener´ atorok tesztel´ese • Brown-mozg´ as
Sz´ am. szim. labor 2015 4. ea.
U(0,1) gener´al´asa
T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF
• C´el: f¨uggetlen, egyenletes eloszl´as´ u ,,v´eletlen” sz´amok gener´al´asa [0, 1]-en • Leggyakrabban eg´eszek gener´alva 0 ´es egy maxim´alis ´ert´ek k¨ oz¨ott • De: a sz´am´ıt´ og´ep determinisztikusan m˝ uk¨odik, nincs v´eletlen, csak ´alv´eletlen
Hivatkoz´ asok
• Line´aris kongruenci´alis algoritmus (csak az el˝ oz˝o sz´amot kell t´arolni): Xn+1 = (aXn + c) %m Rendszer glibc (gcc) ANSI C C99, C11 Visual C++
m 231 231 232 232
a 1103515245 1103515245 1103515245 214013
• Gond: peri´odusok, cs´ıkoz´as...
c 12345 12345 12345 2531011
Eredm´eny b30 . . . b0 b30 . . . b16 b30 . . . b16 b30 . . . b16
Sz´ am. szim. labor 2015 4. ea.
U(0,1) gener´al´asa
T˝ oke Csaba U(0,1)
• A Fibonacci-sorozat mint´ aj´ ara:
GSL
Xn = Xn−j ⊗ Xn−k
Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
ahol (j, k) r¨ ogz´ıtett, ⊗ lehet XOR, ¨ osszead´as, b´armi.
• Speci´ alisan: Tausworthe-gener´ ator (a megel˝oz˝o 250
sz´amot kell t´ arolni):
Xn = Xn−250 XOR Xn−103 • Gond: hogyan inicializ´ aljuk a gener´ atort?
´ Orai feladat: ajlban a double mca(void) f¨ uggv´eny Az rng defs orai.c f´ kieg´esz´ıt´ese, hogy a multiplikat´ıv kongruenci´alis m´ odszert haszn´alja. Legyen m = 231 − 1, a = 16807, c = 0, a kimenet [0, 1]-ben. Az el˝ oz˝ o ´ert´eket a f¨ uggv´enynek kell t´arolnia.
Sz´ am. szim. labor 2015 4. ea.
GSL v´eletlensz´am-gener´ator
T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
• v´ altoz´o: gsl rng *r;
• k¨ ornyezeti v´ altoz´ok olvas´ asa (GSL RNG SEED,
GSL RNG TYPE): gsl rng env setup;
• gener´ ator inicializ´ al´ asa:
r=gsl rng alloc(gsl rng default); • gsl rng default==gsl rng mt19937 Mersenne twister
algoritmus, peri´ odus 219937 − 1 • lehet m´ as gener´atort is haszn´alni • seed megad´ asa: gsl rng set(r,seed); • v´ eletlen sz´am: unsigned long gsl rng get(r); • [0, 1]-en: double gsl rng uniform(r); • felszabad´ıt´ as: gsl rng free(r);
Sz´ am. szim. labor 2015 4. ea.
Seed
T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
• Csak egyszer kell megadni, a program elej´ en
• GSL: gsl rng set(r,seed);
• Tipikus haszn´ alat seed=time(NULL); • http://www.random.org/
• A megism´ etelhet˝ os´eg miatt ´erdemes elt´arolni a haszn´alt
seedet
Sz´ am. szim. labor 2015 4. ea. T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
Adott eloszl´as´u v´eletlen sz´amok • Lehet inverz eloszl´ asf¨ uggv´enyb˝ ol (inverse transform
sampling), ha az egyszer˝ uen sz´ amolhat´ o: −1 FX (x) = P (X < x), ekkor X = FX (U), U ∼ U (0, 1), nem mindig m˝ uk¨ odik (diszkr´et eloszl´ asok, norm´alis eloszl´as) • Eloszl´ as transzform´ aci´ oja: p´eld´ aul norm´alis eloszl´asra Box–M¨ u ller-transzform´ a ci´ o , U, V √∼ U (0, 1) f¨ uggetlenek, √ Z1 = −2 ln U cos (2πV ), Z2 = −2 ln U sin (2πV ), ekkor Z1 , Z2 ∼ N (0, 1) f¨ uggetlenek • Visszautas´ıt´ asos mintav´etelez´es (rejection sampling): f (x) s˝ ur˝ us´egf¨ uggv´eny szerint kell sorsolni, g (x) szerint k¨onny˝ u sorsolni (p´eld´ aul egyenletes eloszl´ as) ´es f (x) ≤ Mg (x), M ≥ 1, ekkor sorsolni kell egy sz´amot g (x) szerint, majd ezt elfogadni f (x) szerint sorsoltnak f (x) osz´ın˝ us´eggel, ami egy U ∼ U (0, 1) Mg (x) ≤ 1 val´ 1 sorsol´as´aval megtehet˝ o, de az elfogad´ as csak az esetek M r´esz´eben t¨ort´enik meg, ´ıgy rossz a hat´ asfok
Sz´ am. szim. labor 2015 4. ea.
Adott eloszl´as´u v´eletlen sz´amok
T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
Egyenletes: U(a, b), a < b ∈ R • s˝ ur˝ us´egf¨ uggv´eny: f (x, x ∈ [a, b]) = 2
(b−a) 2 • E(X ) = a+b 2 , D (X ) = 12 • lehet m´ as eloszl´ asok gener´al´as´ara,
tesztel´es´ere haszn´ alni • gsl ran flat(r,a,b);
Bernoulli: Ber (p), p ∈ [0, 1] • P(X = 1) = p, P(X = 0) = 1 − p • E(X ) = p, D2 (X ) = p (1 − p) • p´ elda: p´enzfeldob´ as • gsl ran bernoulli(r,p);
1 b−a
Sz´ am. szim. labor 2015 4. ea.
Adott eloszl´as´u v´eletlen sz´amok
T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
Binomi´alis: Binom(n, p), n ∈ N, p ∈ [0, 1] • P(X = k, k = 0, . . . , n) = kn p k (1 − p)n−k • E(X ) = np, D2 (X ) = np (1 − p) • p´ elda: sok f¨ uggetlen p´enzfeldob´as • gsl ran binomial(r,p,n); ´ Orai feladat: az rng types orai.c f´ ajlban az int binomial(n,p) f¨ uggv´eny meg´ır´ asa Geometriai: Geom(p), p ∈ [0, 1] k • P(X = k, k ∈ Z+ 0 ) = p(1 − p) 1−p 1−p 2 • E(X ) = p , D (X ) = p 2 • p´ elda: h´ anyszor kapunk ´ır´ast a k¨ovetkez˝o fej el˝ ott • gsl ran geometric(r,p)-1;
Sz´ am. szim. labor 2015 4. ea.
Adott eloszl´as´u v´eletlen sz´amok
T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
Poisson: Poi (λ), λ ∈ R+ • • • •
k
λ −λ P(X = k, k ∈ Z+ 0 ) = k! e 2 E(X ) = λ, D (X ) = λ p´elda: radioakt´ıv boml´ asok sz´ama gsl ran poisson(r,lambda);
Exponenci´alis: Exp(λ), λ ∈ R+ • s˝ ur˝ us´egf¨ uggv´eny: f (x, x ∈ R+ ) = λe−λx • E(X ) = λ1 , D2 (X ) = λ12 • p´ elda: ´elettartam, v´ arakoz´asi id˝o • gsl ran exponential(r,1.0/lambda); ´ Orai feladat: az rng types orai.c f´ ajlban a double exponential(lambda) f¨ uggv´eny meg´ır´ asa az inverzi´os m´ odszerrel, ha FX (x, x ∈ R+ ) = 1 − e−λx
Sz´ am. szim. labor 2015 4. ea.
Adott eloszl´as´u v´eletlen sz´amok
T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
• X , Y , Z ∼ U(−1, 1), Egyenletes g¨omb¨on: X2 + Y 2 + Z2 = 1 • gener´ al´ as: Z ∼ U(−1, 1), ϕ ∼ U(0, 1), √ X = √1 − Z 2 cos (2πϕ), Y = 1 − Z 2 sin (2πϕ) • p´ elda: v´eletlen ir´ any sorsol´asa • gsl ran dir3d(r,&x,&y,&z);
Norm´alis: N(µ, σ 2 ), µ ∈ R, σ ∈ R+ • s˝ ur˝ us´egf¨ uggv´eny: • • • • •
− (x−µ) √ 1 e 2σ 2 2πσ2 D2 (X ) = σ 2
f (x, x ∈ R) =
2
E(X ) = µ, p´elda: f¨ uggetlen k´ıs´erletek ´atlaga gsl ran gaussian(r,s)+m; (Box-M¨ uller), gsl ran gaussian ziggurat(r,s)+m; gsl ran gaussian ratio method(r,s)+m;
Sz´ am. szim. labor 2015 4. ea.
V´eletlensz´am-gener´atorok tesztel´ese
T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
• azt kell tesztelni, hogy a gener´ ator ´ altal sorsolt sz´amok
mennyire felelnek meg f¨ uggetlen, egyenletes eloszl´as´ u v´eletlen sz´amoknak
• ´ altal´aban valamilyen m´ as eloszl´ as gener´alva, annak
eloszl´asa ¨osszehasonl´ıtva az elm´eleti v´ arakoz´assal
• k¨ onyvt´ar (GCC-hez): TestU01 (L’Ecuyer, Simard)[2] • p´ elda: rng types plot geometric
• ´ abr´azol´as: • plot for [i=2:6] ’test.out’ using 1:i with points pointtype 7 linecolor i pointsize 0.5 title "generator " .(i-3) • plot for [i=3:6] ’test.out’ using 1:(column(i)-$2) with points pointtype 7 linecolor i pointsize 0.5 title "generator " .(i-3)
Sz´ am. szim. labor 2015 4. ea. T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
Brown-mozg´as • apr´ o r´eszecsk´ek diff´ uzi´os mozg´ as´ at ´ırja le valamilyen
k¨ozegben, sok kicsi, v´eletlen ir´ any´ u er˝ o hat´as´ara
• matematikai defin´ıci´ o egy dimenzi´ oban: • B0 = 0 majdnem biztosan • Bt majdnem biztosan folytonos • Bt − Bs ∼ N(0, t − s), t ≥ s ≥ 0 • 0 < t1 < · · · < tn eset´ en Bt1 , Bt2 − Bt1 , . . . , Btn − Btn−1 f¨uggetlenek • Bt s˝ ur˝ us´egf¨ uggv´enye
f (x; t) = √
1 − x2 e 2t 2πt
• ez a ∂t p = 12 ∂x2 p diff´ uzi´os (Fokker–Planck-)egyenlet
megold´asa p(x; 0) = δ(x) kezdeti felt´etellel, v´egtelenben elt˝ un˝o megold´ asokra • sztochasztikus differenci´ alegyenletek (p´eld´aul Langevin-egyenlet) megold´ as´ ara haszn´ alj´ak
Sz´ am. szim. labor 2015 4. ea. T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
Brown-mozg´as • Brown-mozg´ as gener´ al´ asa:
Bt =
Z
t
dBs = lim 0
N→∞
t0 =0,t XN =t
i =0,...,N−1
Bti +1 − Bti
• f¨ uggetlen norm´ alis eloszl´ asok ¨ osszege:
Bti +1 − Bti ∼ N(0, ti +1 − ti ) • vizsg´ alhat´o a p´ alya egy megval´ os´ıt´ asra, illetve a v´ arhat´o ´ert´ek ´es a sz´ or´ as az id˝o f¨ uggv´eny´eben • p´ elda: Brownian motion • ´ abr´azol´as: • plot ’test.out’ every :::0::0 with points
pointtype 7 linecolor 3 pointsize 0.5 title "Sample path" • plot ’test.out’ ev :::1::1 u 1:2 w p pt 7 lc 3 ps 0.5 t "Mean", \ ’test.out’ ev :::1::1 u 1:3 w p pt 7 lc 7 ps 0.5 t "Standard deviation"
Sz´ am. szim. labor 2015 4. ea. T˝ oke Csaba U(0,1) GSL Adott eloszl´ as
H´azi feladat • Brown-h´ıd (Brownian bridge) gener´ al´ asa
Xt
= (1 − t)
Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
=
Z
0
t
1 dBs = 1−s
lim (1 − t)
N→∞
t0 =0,t XN =t
i =0,...,N−1
Bti +1 − Bti , t ∈ [0, 1] 1 − ti
• Vizsg´ aljuk a p´ aly´ at, a v´ arhat´o ´ert´eket ´es a sz´or´ast az id˝o
f¨ uggv´eny´eben. • Bek¨ uldend˝ o a forr´ ask´od (Brownian bridge) ´es h´arom ´abra: egy a p´ aly´ ar´ ol (sample.jpg), egy a v´arhat´o ´ert´ekr˝ol (mean.jpg), egy a sz´ or´ asr´ ol (stddev.jpg), mindegyik a [0, 1] intervallumon az id˝o f¨ uggv´eny´eben, megfelel˝o felbont´assal, az ´ atlagolt mennyis´egekre megfelel˝o sz´am´ u (≈ 1000) ´atlagol´ assal. • Maxim´ alis pontsz´ am: 2 pont
Sz´ am. szim. labor 2015 4. ea.
Hivatkoz´asok
T˝ oke Csaba U(0,1) GSL Adott eloszl´ as Tesztel´ es Brown-mozg´ as HF Hivatkoz´ asok
[1] Gnu Scientific Library Reference Manual. http://www.gnu.org/software/gsl/manual/html node/ [2] Pierre L’Ecuyer and Richard Simard. TestU01 A Software Library in ANSI C for Empirical Testing of Random Number Generators. http://www.iro.umontreal.ca/˜simardr/testu01/tu01.html ´ [3] Ertelmes honlap a Brown-h´ıd matematik´aj´ahoz: www.math.uah.edu/stat/brown/Bridge.html
