Całkowanie metodą Monte Carlo w 2D z rozkładem eksponencjalnym.
Tomasz Chwiej 16 czerwca 2021
1 Opis problemu
Należy oszacować wartość całki
I =
∫1 0
∫1 0
g(x, y)· f(x) · f(y)dxdy (1)
gdzie:
g(x, y) = sin(x + y)
ln(2 + x + y), f (x) = e−x, f (y) = e−y (2) metodą Monte Carlo. Jako dokładną wartość przyjmiemy Idok = 0.2557840245.
1.1 Sposób 1
Używamy bezpośrednio kwadratury Monte Carlo dla funkcji podcałkowej
I ≈ 1 n
∑n i=1
z(xi, yi) = 1 n
∑n i=1
g(xi, yi)· f(xi)· f(yi) (3)
gdzie xi, yi ∈ U(0, 1), a n jest liczbą losowań.
W celu kontroli jakości rozwiązania obliczamy wariancję pojedynczego wyniku:
σ2 = 1 n− 1
∑n
i=1
z2(xi, yi)− 1 n
( n
∑
i=1
z(xi, yi) )2
(4)
oraz odchylenie standardowe wartości średniej:
σsr= σ
√n (5)
1.2 Sposób 2
Zauważmy, że zmienne x oraz y mogą być opisane funkcjami gęstości prawdopodobieństwa f (x) i f (y) (jeszcze nieunormowanymi) które występują w całkowaniu
I =
∫1 0
∫1 0
g(x, y)
CxCy · Cxf (x)· Cyf (y)dxdy (6) Zatem możemy do wyznaczenia całki możemy użyć także wzoru:
I ≈ 1 n
∑n i=1
g(xi, yi) 1
CxCy (7)
1
gdzie xi, yi ∈ Q, gdzie Q jest rozkładem eksponencjalnym, n jest liczbą losowań, a stałe Cx i Cy wyznaczymy dla rozkładu Q. Wariancję pojedynczej zmiennej losowej g(x, y) wyznaczamy ze wzoru
σ2= 1 n− 1
∑n
i=1
(g(xi, yi) CxCy
)2
− 1 n
( n
∑
i=1
g(xi, yi) CxCy
)2
(8)
Generator dla zmiennej o rozkładzie eksponencjalnym skonstruujemy wykorzystując metodę odwra- cania dystrubuanty. Najpierw zajmiemy się stałą normalizacji rozkładu:
fgp(x) = Cxe−x →
∫ 1
0
fgp(x)dx = 1 (9)
skąd otrzymujemy
Cx = Cy = 1
1− e−1 (10)
Zgodnie z metodą odwracania dystrybuanty mamy warunek:
F (x) =
∫ x
0
fgp(x0)dx0 = Cx(1− e−x) = XU ∈ U(0, 1) (11) który pozwala wyznaczyć zmienną XQ o rozkładzie eksponencjalnym Q:
XQ=−ln (
1−XU
Cx
)
(12)
2 Zadania do wykonania
1. Przyjmujemy N = 105 jako maksymalną liczbę losowań.
2. Oszacować wartość całki (1) wykorzystując wzór (3) oraz wyznaczyć błąd oszacowania wartości średniej (wykorzystując wzór 5). W trakcie obliczeń do pliku zapisać: n (aktualna ilość wykona- nych losowań), I(n) oraz σsr(n) dla n = 10, 102, 103, 104, 105.
3. Oszacować wartość całki (1) wykorzystując wzór (7) oraz wyznaczyć błąd oszacowania wartości średniej (wykorzystując wzór 8). W trakcie obliczeń do pliku zapisać: n (aktualna ilość wykona- nych losowań), I(n) oraz σsr(n) dla n = 10, 102, 103, 104, 105.
4. Testowanie generatora Q.
a) Przedział [0, 1] podzielić na k = 10 podprzedziałów, następnie wylosować N = 105 liczb pseudolosowych tym generatorem.
b) Dla każdego podprzedziału wyznaczyć ilość liczb ni, i = 1, 2, . . . , k, które wpadają do podprzedziału i.
c) Wyznaczyć teoretyczne prawdopodobieństwo wylosowania liczby z danego popdprzedziału jako
Pi(xi−1< x¬ xi) = F (xi)− F (xi−1) (13) gdzie dystrubuanta F (x) jest określona wzorem (11), a xi−1 oraz xi wyznaczają krańce i− tego podprzedziału.
d) Narysować histogramy pi oraz Pi na jednym rysunku, np. w Gnuplocie plot ’hist_pk.dat’ u 1:2 w boxes lc rgb ’black’ t ’num’,\
’hist_Pk.dat’ u 1:3 w boxes lc rgb ’red’ t ’dok’
e) Obliczyć wartość statystyki χ2k−1:
χ2k−1 =
∑k i=1
(ni− N · Pi)2 N · Pi
(14) i zapisać do pliku.
3 Przykładowe wyniki
0.252 0.254 0.256 0.258 0.26 0.262 0.264 0.266 0.268 0.27 0.272 0.274
101 102 103 104 105
Inum
N
sposob-1 sposob-2
10-5 10-4 10-3 10-2 10-1
101 102 103 104 105
error
N
σ1
|Idok-Inum|1 σ2
|Idok-Inum|2
Rysunek 1: Wartość całki (lewy) i błąd (prawy) liczona 1 i 2 sposobem.
0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15 0.16
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 pi
x
num dok
Rysunek 2: Histogram dla generatora eksponencjalnego Q.
