Metoda Monte Carlo – cvičení

Lineární kongruenční generátory

Nechť p je prvočíslo, x0p1^ a a je generátor multiplikativní grupy modulo p. Definujme rekurentní posloupnost xnaxn1modp. Potom tato posloupnost bude postupně procházet všechna čísla z p1^ v pseudonáhodném pořadí. Uděláme-li z ní náhodnou proměnnou X, potom bude XU(p1^) neboli pro dostatečně velké p máme YXpU((0,1)). Ovšem zkusíme-li vyplnit čtverec tím, že střídavě generujeme souřadnice x a y, zjistíme, že ho to nevyplní, ale bude to na něm vytvářet jakési čáry.

Spočtěme si pro Y střední hodnotu a rozptyl:

EY=01xfY(x)dx=01xdx=12,
VarY=E(X12)=112.

Ze zákona velkých čísel víme, že definujeme-li Zi=1mYi pro dostatečně velké m, bude přibližně ZN(m2,m12). Když z toho chceme udělat standardní normální rozdělení, znormujeme to:

QZm2m12.

Ale pozor, zvolíme-li pro jednoduchost třeba m=12, tahle náhodná proměnná ve skutečnosti vždy bude v intervalu (6,6). Tedy nemáme skutečně normální rozdělení.

Výpočet objemu

Pojďme si metodou Monte Carlo spočítat objem paraboloidu. Mějme bod u=(x,y,z)[R,R]2×[0,H] a paraboloid A daný nerovnicí zHR2(x2+y2). Budeme-li bod brát rovnoměrně náhodně, potom objem paraboloidu můžeme spočítat jako

V=4R2Hp[uA]p.

Samozřejmě abychom tuto pravděpodobnost spočetli, musíme hodněkrát náhodně střílet. Budeme tedy mít náhodnou proměnnou MBi(n,p). Poté můžeme odhadnout pMN. Skutečně platí EM=Np. Ale s jakým rozptylem?

σ2=Np(1p)NMN(1MN)=M(1MN)
σ=M1MN
γ=σEM=1MNM<1M=𝒪︀(M12)=𝒪︀(N12)

Taky můžeme použít polární souřadnice. Zvolíme náhodně φU(0,2π). Ovšem u poloměru už nechceme rovnoměrné rozdělení, protože bodů vzdálenějších od středu je víc. Na přednášce si odvodíme, že je vhodné zvolit rU(0,1). Potom máme nerovnici zHR2r a stačí spočítat πR2Hp.

Ještě jiný přístup. Místo abychom náhodně stříleli body v trojrozměrném prostoru, je budeme střílet ve dvojrozměrném a pro každý spočteme výšku paraboloidu v daném bodě. Tuto výšku spočteme jako

z=(1r2R2)H=(1x2+y2R2)Hg(x,y),
Eg(X,Y)1Ni=1Ng(xi,yi)g¯,
sx2=1N1i=1N(g(xi,yi)g¯)2.

Z toho si už můžeme snadno spočíst objem:

V=πr2Eg(X,Y).

Nerovnoměrná rozdělení

Chceme generovat náhodnou proměnnou s exponenciální hustotou pravděpodobnosti
f(x)=[x0]1μexp(xμ).
Odpovídající distribuční funkce je
F(x)=xf=[x0](1exp(xμ)).
Jak takovéto rozdělení nagenerujeme z náhodné proměnné s rovnoměrným rozdělením rR(0,1)? Stačí vyřešit rovnici F(x)=r. Vyjde nám
x=μln(1r)=μlnr*,r*R(0,1).
Tento přístup můžeme samozřejmě použít i pro jiná rozdělení. Mějme třeba Cauchyovo rozdělení
f(x)=1π(1+x2),
F(x)=12+arctanxπ.
Opět budeme řešit rovnici F(x)=r, čímž dostaneme
x=tan(π2(2r1)).
Jako poslední příklad si zkusme tímto způsobem vygenerovat rovnoměrné rozdělení xR(a,b).
f(x)=[x[a,b]]ba,
F(x)={0x<a,xabax[a,b],1x>b.
Řešením rovnice F(x)=r, přičemž uvažujeme jenom „prostou část“ uprostřed, dostaneme
x=a+(ba)r.
To jsme samozřejmě mohli uhodnout, ale je dobré si ověřit, že systematický přístup sedí s intuicí.

Testování náhodné volby

Mějme nějaké náhodné rozdělení. Zvolíme si body x1<x2<<xn1, přičemž bereme x0,xn. Tím si rozdělíme reálnou osu na šuplíčky. Pravděpodobnost, že se náhodná proměnná bude vyskytovat v k-tém šuplíčku, je
pkxk1xkf=F(xk)F(xk+1).
Provedeme-li N pozorování, dostaneme nějaká čísla N1,,Nn vyjadřující, kolikrát jsme se trefili do kterého šuplíčku. Střední hodnota bude
Ok=Npk.
V praxi chceme šuplíčky zvolit tak, aby bylo Ok5. Podle nějaké Pearsonovy věty platí
χ2k=1n(NkOk)2Okχn12.

Box-Müllerova transformace

Cílem je vygenerovat z rovnoměrného rozdělení normální rozdělení. Kdybychom chtěli použít obecnou metodu pro nerovnoměrná rozdělení, dopadli bychom blbě, protože jak víme, kumulativní distribuční funkce nemá elementární tvar:
F(x)=xf=x12πexp(t22)dt.
Použijeme oblíbený trik, že si to rozšíříme do dvou rozměrů. Mějme funkci
f2(x,y)12πexp(x2+y22).
Přejdeme-li do polárních souřadnic, máme φU((o,2π)). Pro distribuční funkci poloměru platí
FR(r)=x2+y2r212πexp(x2+y22)d(x,y)=0r02π12πexp(r~22)dφdr~==1exp(r22).
Tuhle funkci už dokážeme invertovat: máme-li rovnoměrnou proměnnou u, potom FR(r)=u, tedy r=2ln(1u)=2lnu*.

Poissonovo rozdělení

Máme exponenciální rozdělení XEx(μ),
f(x)=1μexp(xμ)I(x0),
F(x)=(1exp(xμ))I(x0).
Dále existuje binomické rozdělení YBi(N,p),
py(N,p)=(Ny)py(1p)Ny.
Jeho limitou pro N,Np=λ je Poissonovo rozdělení XPo(λ),
px(λ)=limNpy(N,λN)=λxx!exp(λ).
Teď nás zajímá, jak Poissonovo rozdělení generovat. Mohli bychom opět pro náhodnou proměnnou UU(0,1) řešit rovnici F(x)=u. Nyní už ale F není prostá funkce, takže si musíme ze všech možných řešení zvolit to nejmenší. Pro zjednodušení můžeme počítat
F0=p0=exp(λ),Fx=Fx1+px=Fx1+λxpx1.
Další možnost je aproximovat ho pomocí binomického, kde můžeme počítat
F0=p0=(1p)N,Fy=Fy1+py=Fy1+p(Ny+1)(1p)ypy1.

Výpočet integrálu

Chceme spočítat Iabf(x)dx. Víme, li, že Zf(x)Q, můžeme vzít náhodné proměnné XU(a,b),YU(Z,Q) a máme
P[Yf(X)]=ab(f(x)z)dx(ba)(Qz)pHIT,
I=(ba)(Z+(QZ)pHIT).
Střelíme-li náhodně N-krát a z toho se K-krát trefíme, můžeme aproximovat pHITKN. Pro tuto metodu je
σ=(QZ)pHIT(1pHIT)N=𝒪︀(1N)
Jiný způsob, jak by na to šlo jít, je vzít si náhodnou proměnnou XU(a,b) a počítat
Ef(X)=1baabf(x)dx=Iba,
I=(ba)Ef(X)baNi=1Nf(Xi),xiU(a,b).
Pro tuto metodu je
σ=1N(ba)ab(g(x))2dx(abg(x)dx)2=𝒪︀(N).
Rozptyl je menší než u předchozí metody.

Nevlastní integrál

Zkusme spočítat integrál
0cosωx1+x2dx.
Potřebujeme zavést nějakou substituci, která to dostane do omezeného intervalu. Jedna možnost je ξx1+x. Tím dostaneme
01cosωξ1ξ(1ξ)2+ξ2dξ.
Další možnost by byla ξexp(x), čímž dostaneme
01cos(ωlnξ)(1+ln2ξ)ξdξ.
Nebo taky ξ2πarctanx:
01πcos(ωtanπξ2)(1+tan2πξ2)2cos2πξ2.

Vícerozměrné integrály

Ukážeme si, v čem metoda Monte Carlo vyniká oproti deterministickým metodám integrace. Chceme-li spočítat abf dělením intervalu na N kousíčků velikosti habN a chceme přesnost ε=𝒪︀(hk), kde k je řád konvergence metody, časová náročnost bude
T(N)=𝒪︀(N)=𝒪︀(h1)=𝒪︀(ε1k).
Počítáme-li d-rozměrný integrál, kde j-tou dimenzi dělíme na Nj kousíčků, časová náročnost bude
T(N)=𝒪︀(j=1dNj)=𝒪︀(j=1dhj1)=𝒪︀(hd)=𝒪︀(εdk).
Kdybychom to na oblasti Od počítali metodou Monte Carlo, tak jak už jsme si odvozovali, přesnost bude
ε=σV(O)T(N),
z čehož máme T(N)=𝒪︀(ε2). Metoda Monte Carlo bude tedy rychlejší v případě, kdy dk>2 neboli d>2k.