Mamy dany wektor losowy \((X, Y)\) o standardowym dwuwymiarowym rozkładzie normalnym, gdzie \(X\) i \(Y\) są niezależne. Stąd łączna gęstość jest dana wzorem \[f_{X,Y}(x,y)=\frac{1}{2\pi}e^{-\frac12(x^2+y^2)}, \ x,y\in \mathbb{R}.\] Ponieważ zmienne losowe są niezależne, to rozkład łączny jest iloczynem rozkładów brzegowych, tj. \(f_{X,Y}(x,y)=f_X(x)f_Y(y)\). Zadanie sprowadza się do obliczenia odpowiednich prawdopodobieństw i aby je wyznaczyć skorzystamy z twierdzenia

Twierdzenie 1

Jeśli \(X,Y\) są ciągłymi zmiennymi losowymi, to prawdopodobieństwo, że \((X,Y) \in D, \ D\subseteq \mathbb{R}^2\) wyraża się wzorem \[P((X,Y) \in D)=\int\!\!\!\int_{D} f_{X,Y}(x,y)dxdy.\]

Następnie policzymy odpowiednie całki, korzystając z metod numerycznych, a mianowicie aproksymacji Riemanna i metody Monte Carlo. Na końcu dokonamy porównania metod.

Przybliżmy po krótce obie metody.

METODA APROKSYMACJI RIEMANNA

Metoda Riemanna w \(R^1\) polega na przybliżaniu pola pod wykresem polami prostokątów. W analogicznej wersji, w przypadku \(R^2\), objętość między funkcją a obszarem całkowanie, przybliżamy licząc objętości prostopadłościanów. Dla uproszczenia rozważmy obszar całkowania, który jest kwadratem. Mamy policzyć całkę \(\int_a^b\int_{a}^b g(x,y)dxdy\).

Zatem zaczynamy od podziału kwadratu \([a,b]\times [a,b]\) na mniejsze kwadraty, dzieląc \(x\in[a,b]\) na \(n\) odcinków i analogicznie \(y\). Stąd dostajemy \(n^2\) kwadratów o wymiarach \(\frac{(b-a)}{n}\times \frac{(b-a)}{n}\). Należy policzyć wartość funkcji nad każdym kwadratem, w punktach środkowych kwadratu \(x_{s_i}=\frac{x_{i+1}-x_{i}}{2}\), gdzie \(i=0,1,2,...,n\) i \(y_{s_i}=\frac{y_{i+1}-y_{i}}{2}\), gdzie \(i=0,1,2,...,n\). Możemy przybliżyć całkę, korzystając ze wzoru \[\int_a^b\int_{a}^b g(x,y)dxdy \approx \frac{(b-a)}{n} \frac{(b-a)}{n}\sum_{i=0}^n\sum_{j=0}^n g(x_{s_{i}},y_{s_{j}}).\] Przy czym przybliżenie jest tym większe, im większe jest \(n\).

Gdy obszar nie jest prostokątem, ale jest obszarem normalnym, dokonujemy podziału kwadratu zawierającego ten obszar a następnie wybieramy punkty spełniające równanie obszaru.

METODA MONTE CARLO

W przeciwieństwie do aproksymacji Riemanna metoda Monte Carlo korzysta z metod rachunku prawdopodobieństwa. Skupimy się na przedstawieniu metody Monte Carlo do całkowania funkcji dwóch zmiennych. Załóżmy, że mamy policzyć całkę \(\int_a^b\int_{a}^b g(x,y)dxdy\), znów dla uproszczenia opisu na kwadracie.

Jeśli mamy wektor losowy \((X,Y)\) o gęstości łącznej \(f(x,y)\), to wartość oczekiwaną możemy policzyć z definicji \[E[g(X,Y)]=\int\!\!\!\int_{\mathbb{R}}g(x,y)f(x,y)dxdy.\] Rozważmy próby losowe \(X_1,...,X_n\) oraz \(Y_1,...,Y_n\) o rozkładach jednostajnych na \([a,b]\). Gdy zmienne losowe \(\xi_i,\eta_i\in U(0,1)\) dla \(i=1,...,n\), to \[X_i=a+\xi_i(b-a), \ \ \ \ Y_i=a+\eta_i(b-a).\] Ponieważ zmienne losowe są niezależne, to ich łączna gęstość jest iloczynem gęstości brzegowych i \(f(x,y)=(\frac{1}{b-a})^2\). Stąd \[E[g(X,Y)]=\int\!\!\!\int_{\mathbb{R}}g(x,y)f(x,y)dxdy\approx \frac{(b-a)^2}{n}\sum_{i=1}^ng(X_i,Y_i).\] Z Mocnego Prawa Wielkich Liczb wiemy, że \[\lim_{n\to\infty}\frac{(b-a)^2}{n}\sum_{i=1}^ng(X_i,Y_i)=E[g(X,Y)] \ \ \text{ z pr. 1}\] Zatem zbieżność metody zapewnia nam MPWL. Dodatkowo przybliżenie jest dokładniejsze, im większe jest \(n\). Nad tempem zbieżności i poprawą dokładności zastanowimy się pod koniec.

Przedstawimy teraz instrukcję postępowania do obliczenia całki podwójnej na obszarze \(D\) metodą Monte Carlo:

Mając za sobą wstęp teoretyczny, możemy przejść do obliczenia zadanych prawdopodobieństw. Całki zostały policzone numerycznie w programie R. Każde z trzech prawdopodobieństw wyznaczymy trzema sposobami: analitycznie, aproksymacją Riemanna i metodą Monte Carlo. Następnie porównamy zbieżność i dokładność obu numerycznych metod.

Wyliczenie \(P(0< X\leq1, 0< Y\leq1)\)

Z Twierdzenia 1 dane prawdopodobieństwo możemy policzyć jako \[P(0< X\leq1, 0< Y\leq1)= \int_0^1\int_0^1 \frac{1}{2\pi}e^{-\frac12(x^2+y^2)} dxdy.\] Powyższej całki nie da się policzyć w sposób analityczny, dlatego skorzystajmy z niezależności \(X\) i \(Y\), i policzmy to prawdopodobieństwo jako iloczyn prawdopodobieństw i odczytamy wartości dystrybuanty rozkładu normalnego z tablic.

\[P(0< X\leq 1, 0< Y\leq1)=P(0< X\leq1)P(0< Y\leq 1)=[\phi(1)-\phi(0)]^2 \approx 0,116516.\]

Przyjrzyjmy się wynikom obliczeń całki w sposób numeryczny. Najpierw metoda aproksymacji Riemanna.

Zgodnie z opisem teoretycznym dokonaliśmy podziału obszaru całkowania, kwadratu, w naszym przypadku na 10000 kwadratów (oś X i Y dzielimy na m=100 odcinków). Następnie obliczamy wartość funkcji dla każdego kwadratu, ponieważ podział jest dość gęsty nie liczymy wartości w środku kwadratu, a w wierzchołkach. Nie ma to dużego wpływu na wartość całki.

#a) calka na kwadracie metoda aproksymacji Riemanna
f <- function(x,y) 1/(2*pi)*exp(-0.5*(x^2+y^2))  #funkcja gestosci dwuwymiarowygo rozkladu normalnego gdy \rho=0
m=100 #gestosc podzialu obszaru
a=0
b=1
#podzial obszaru na kwadraty
x=seq(a,b,length.out=m)
y=seq(a,b,length.out=m)
punkty=expand.grid(x,y) #macierz wspólrzednych kwadratow
#obliczanie calki podwójnej: pole obszaru*wartosci funkcji w wylosowanych punktach 
p=((b-a)/m)^2
sum(p*f(punkty[1],punkty[2]))
## [1] 0.1163723

W wyniku tego całkowania otrzymujemy wynik \[ \int_0^1\int_0^1 \frac{1}{2\pi}e^{-\frac12(x^2+y^2)} dxdy \approx 0,1163723.\] Zobaczmy teraz wynik metody Monte Carlo

#a) calka na kwadracie Monte Carlo
f <- function(x,y) 1/(2*pi)*exp(-0.5*(x^2+y^2))  #funkcja gestosci dwuwymiarowygo rozkladu normalnego gdy \rho=0
m = 10000 #liczba wylosowanych punktów

a=0
b=1
#losowanie liczb z U(a,b)
x= a+(b-a)*runif(m)
y= a+(b-a)*runif(m)
plot(x,y) #wykres wylosowanych punktów z obszaru calkowania

#obliczanie calki podwójnej: pole obszaru*wartosci funkcji w wylosowanych punktach 
(b-a)^2*mean(f(x,y))
## [1] 0.1165156

W wyniku tego całkowania otrzymujemy wynik \[ \int_0^1\int_0^1 \frac{1}{2\pi}e^{-\frac12(x^2+y^2)} dxdy \approx 0,1161215.\] Jak widzimy obie metody dają taki sam rezultat do trzeciego miejsca po przecinku. Jest to dość dobre przybliżenie, jednak należy pamiętać, że za każdym puszczeniem program, w metodzie Monte Carlo, losuje nowe liczby, przez co wartości całki wahają się nieznacznie po trzecim miejscu po przecinku.

Porównajmy błędy względne obu metod, mając na uwadze wahania w metodzie MC.

#porównanie obu metod ze wzgledu na liczbe m 
a=0
b=1
c1={}
c2={}
n <- c(100,1000,10000,100000,1000000) #liczba losowan punktu\liczba kwadratow w podziale
for (i in 1:5) {
  # MC
  x= a+(b-a)*runif(n[i])
  y= a+(b-a)*runif(n[i])
  c1[i] = (b-a)^2*mean(f(x,y))  #wektor calki MC dla różnych n
  # R
  x=seq(a,b,length.out=floor(sqrt(n[i]))) #bierzemy pierwiastek z n, bo chcemy mieć n kwadratów
  y=seq(a,b,length.out=floor(sqrt(n[i])))
  punkty={}
  punkty=expand.grid(x,y)
  p=((b-a)/floor(sqrt(n[i])))^2
  c2[i]=sum(p*f(punkty[1],punkty[2])) #wektor calki R dla różnych n
  
}

wartosc = 0.116516 # najlepsze przybliżenie wartosci calki
blad1 = abs(c1-wartosc)/wartosc #blad wzgledny MC
blad2 = abs(c2-wartosc)/wartosc #blad wzgledny R
plot(c(2,3,4,5,6), blad1, type = 'o', col = 'red',ylim = c(0,0.05), 
     xlab = '10^n podzialu odc/losowań', ylab = 'blad wzgledny przybliżenia', main = 'Porównanie bledów metody R i MC na kwadracie')
  lines(c(2,3,4,5,6),blad2,type = 'o', col= 'blue')

Czerwona krzywa - wykres błędów względnych MC, niebieska - błędy R. Liczba losowań/podziału, to \(10^2,10^3,10^4,10^5,10^6\).

Z wykresu możemy odczytać, że przy podziale/losowaniu dla \(m=1000\) błąd względny jest na poziomie 0,5%. Na kwadracie nie ma dużych różnic między metodami.

Wyliczenie \(P(X^2 + Y^2<1)\)

Zacznijmy od wyliczeń analitycznych \[P(X^2 + Y^2<1)=\int\int_{K(0,1)} \frac{1}{2\pi}e^{-\frac12(x^2+y^2)} dxdy=\] Zastosujemy współrzędne biegunowe

\(x=r\cos\phi, \ y=r\sin\phi, \ x^2+y^2=r^2, \ |J|=r.\) Granice: \(r\in (0,1], \ \phi \in(-\pi,\pi]\). \[=\int_0^1\int_{-\pi}^{\pi} \frac{1}{2\pi}re^{-\frac12 r^2} drd\phi=\int_0^1r e^{-\frac12 r^2}dr=-e^{-\frac12 r^2}\Big|_0^1=1-\frac{1}{\sqrt e}\approx 0,3934693.\]

Możemy przeanalizować teraz metodę Riemanna

#b) calka po kole metoda aproksymacji Riemanna
m=100 #gestosc podzialu obszaru
a=-1
b=1
#podzial obszaru na kwadraty
x=seq(a,b,length.out=m)
y=seq(a,b,length.out=m)
punkty=expand.grid(x,y)
kolo=punkty[punkty[1]^2+punkty[2]^2<1,] #obszar calkowania
plot(kolo,xlim=c(-1,1),ylim=c(-1,1)) #wykres obszaru

#pole jednego kwadratu
p=(b-a)^2/(m^2)
#calka
sum(p*f(kolo[1],kolo[2]))
## [1] 0.3844913

Jedyną różnicą w kodzie w stosunku do podpunktu a) jest obszar całkowania. Musimy brać te kwadraty (punkty wyznaczające kwadraty), które należą do koła o środku 0 i promieniu 1. Otrzymujemy wynik \[\int\int_{K(0,1)} \frac{1}{2\pi}e^{-\frac12(x^2+y^2)} dxdy \approx 0,3844913.\] Teraz metodą Monte Carlo

#b) calka po kole Monte Carlo
f <- function(x,y) 1/(2*pi)*exp(-0.5*(x^2+y^2))  #funkcja gestosci dwuwymiarowygo rozkladu normalnego gdy \rho=0
m = 10000 #liczba wylosowanych punktów

a=-1
b=1
x= a+(b-a)*runif(m)
y= a+(b-a)*runif(m)
#pole obszaru calkowania
Pole_ob=(b-a)^2*length(x[x^2+y^2<1])/length(x)
plot(x[x^2+y^2<1],y[x^2+y^2<1])  #wykres wylosowanych punktów z obszaru calkowania

#calka
Pole_ob*mean(f(x[x^2+y^2<1],y[x^2+y^2<1]))
## [1] 0.3930216

i otrzymujemy wynik \[\int\int_{K(0,1)} \frac{1}{2\pi}e^{-\frac12(x^2+y^2)} dxdy \approx 0,3920665.\] Możemy zauważyć, że wartość całki wyznaczona metodą Monte Carlo jest dokładniejsza.

Zobaczymy porównanie błędów względnych dla obu metod ze względu na \(m\).

#porównanie obu metod ze wzgledu na liczbe m 
a=-1
b=1
c1={}
c2={}
n <- c(100,1000,10000,100000,1000000) #liczba losowan punktu\liczba kwadratow w podziale
for (i in 1:5) {
  # MC
  x= a+(b-a)*runif(n[i])
  y= a+(b-a)*runif(n[i])
  Pole_ob=(b-a)^2*length(x[x^2+y^2<1])/length(x)
  c1[i] =Pole_ob*mean(f(x[x^2+y^2<1],y[x^2+y^2<1]))  #wektor calki MC dla różnych n
  # R
  x=seq(a,b,length.out=floor(sqrt(n[i]))) #bierzemy pierwiastek z n, bo chcemy mieć n kwadratów
  y=seq(a,b,length.out=floor(sqrt(n[i])))
  punkty={}
  punkty=expand.grid(x,y)
  kolo=punkty[punkty[1]^2+punkty[2]^2<1,]
  p=((b-a)/floor(sqrt(n[i])))^2
  c2[i]=sum(p*f(kolo[1],kolo[2])) #wektor calki R dla różnych n
  
}

wartosc = 1-sqrt(exp(-1)) # najlepsze przybliżenie wartosci calki
blad1 = abs(c1-wartosc)/wartosc #blad wzgledny MC
blad2 = abs(c2-wartosc)/wartosc #blad wzgledny R
plot(c(2,3,4,5,6), blad1, type = 'o', col = 'red',ylim = c(0,0.25), 
     xlab = '10^n podzialu odc/losowań', ylab = 'blad wzgledny przybliżenia', main = 'Porównanie bledów metody R i MC na kole')
lines(c(2,3,4,5,6),blad2,type = 'o', col= 'blue')

Czerwona krzywa - wykres błędów względnych MC, niebieska - błędy R. Liczba losowań/podziału, to \(10^2,10^3,10^4,10^5,10^6.\)

Widzimy, że w tym wypadku metoda MC wypada lepiej - generuje mniejszy błąd. Można to wytłumaczyć tym, że obszaru całkowania - koła - nie można podzielić na kwadraty. Dopiero przy dużym podziale, \(m\geq 100000\), kwadraty z podziału “nie wystają” poza obszar.

Wyliczenie \(P(X>0, Y>0, X + Y<1)\)

Całki nie da się wyliczyć analitycznie, dlatego skorzystamy z przybliżenia programu WolframAlpha: \[P(X>0, Y>0, X + Y<1)=\int^1_0\int_{0}^{1-x} \frac{1}{2\pi}e^{-\frac12(x^2+y^2)} dxdy \approx 0,06773. \] Kod programu:

#c) calka po trójkacie metoda aproksymacji Riemanna
m=100 #gestosc podzialu obszaru
a=0
b=1
#podzial obszaru na kwadraty
x=seq(a,b,length.out=m)
y=seq(a,b,length.out=m)
punkty=expand.grid(x,y)
trojkat=punkty[punkty[1]+punkty[2]<1,] #obszar calkowania
plot(trojkat,xlim=c(0,1),ylim=c(0,1)) #wykres obszaru

#pole jednego kwadratu
p=(b-a)^2/(m^2)
#calka
sum(p*f(trojkat[1],trojkat[2]))
## [1] 0.0671654

Aproksymacja Riemanna \[\int^1_0\int_{0}^{1-x} \frac{1}{2\pi}e^{-\frac12(x^2+y^2)} dxdy \approx 0,0671654.\]

#c) calka po trójkacie Monte Carlo
f <- function(x,y) 1/(2*pi)*exp(-0.5*(x^2+y^2))  #funkcja gestosci dwuwymiarowygo rozkladu normalnego gdy \rho=0
m = 10000 #liczba wylosowanych punktów

a=0
b=1
#losowanie liczb z U(a,b)
x= a+(b-a)*runif(m)
y= a+(b-a)*runif(m)
#pole obszaru calkowania
Pole_ob=(length(x[y+x<1]))/length(x)
plot(x[y+x<1],y[y+x<1]) #wykres wylosowanych punktów z obszaru calkowania

#calka
Pole_ob*mean(f(x[y+x<1],y[y+x<1]))
## [1] 0.06749742

Metoda Monte Carlo \[\int^1_0\int_{0}^{1-x} \frac{1}{2\pi}e^{-\frac12(x^2+y^2)} dxdy \approx 0,0678756.\] Przedstawmy od razu wykres błędów

#porównanie obu metod ze wzgledu na liczbe m 
a=0
b=1
c1={}
c2={}
n <- c(100,1000,10000,100000,1000000) #liczba losowan punktu\liczba kwadratow w podziale
for (i in 1:5) {
  # MC
  x= a+(b-a)*runif(n[i])
  y= a+(b-a)*runif(n[i])
  Pole_ob=(length(x[y+x<1]))/length(x)
  c1[i] =Pole_ob*mean(f(x[x+y<1],y[x+y<1]))  #wektor calki MC dla różnych n
  # R
  x=seq(a,b,length.out=floor(sqrt(n[i]))) #bierzemy pierwiastek z n, bo chcemy mieć n kwadratów
  y=seq(a,b,length.out=floor(sqrt(n[i])))
  punkty={}
  punkty=expand.grid(x,y)
  trojkat=punkty[punkty[1]+punkty[2]<1,]
  p=((b-a)/floor(sqrt(n[i])))^2
  c2[i]=sum(p*f(trojkat[1],trojkat[2])) #wektor calki R dla różnych n
  
}

wartosc =0.06773 # najlepsze przybliżenie wartosci calki
blad1 = abs(c1-wartosc)/wartosc #blad wzgledny MC
blad2 = abs(c2-wartosc)/wartosc #blad wzgledny R
plot(c(2,3,4,5,6), blad1, type = 'o', col = 'red',ylim = c(0,0.25), 
     xlab = '10^n podzialu odc/losowań', ylab = 'blad wzgledny przybliżenia', main = 'Porównanie bledów metody R i MC na trójkacie')
lines(c(2,3,4,5,6),blad2,type = 'o', col= 'blue')

Czerwona krzywa - wykres błędów względnych MC, niebieska - błędy R. Liczba losowań/podziału, to \(10^2,10^3,10^4,10^5,10^6\).

Widzimy, że tu, podobnie jak w podpunkcie a, całki wyliczone obiema metodami dają wyniki o podobnej dokładności.

WNIOSKI

Metoda aproksymacji Riemanna i metoda Monte Carlo dają podobne rezultaty na regularnych, prostokątnych obszarach całkowania. MC daje różne wyniki przy każdorazowym wyliczaniu całki, przez co dla małych \(m\) może czasem generować większy błąd. Z tego powodu bezpieczniejszą metodą jest metoda Riemanna. Gdy obszar całkowania nie jest regularny lub nie jest ograniczony prostymi, MC wypada znacznie lepiej.

POPRAWA ZBIEŻNOŚCI ALGORYTMU MC

Możemy teraz zadać sobie pytanie, czy jesteśmy w stanie uzyskać dokładniejszy wynik. Z treści MPWL wynika, że przybliżenie zbiega do faktycznej wartości przy \(m\to\infty\), zatem naturalnym jest próbowanie zwiększenia dokładności poprzez zwiększanie ilości punktów tj. \(m\). Jednakże z Centralnego Twierdzenia Granicznego wiemy, że odchylenie standardowe estymatora jest proporcjonalne do \(\frac{1}{\sqrt{m}}\), zatem \(\sigma(\cdot)=\frac{c}{\sqrt{m}}\), gdzie \(c\) jest stałą. Z tego wynika, że aby zmniejszyć błąd o połowe należy wziąć poczwórnie dużą próbkę. Istnieją inne metody redukcji wariancji w celu poprawienia dokładności i my przyjrzymy się jednej, najprostszej metodzie - metodzie zmiennych antytetycznych.

Mamy policzyć całkę \(\int\!\!\!\int_{D}g(x,y)dxdy\) i wiemy, że sprowadza się to do policzenia \(E|g(X,Y)|\). Przypuśćmy, że mamy pary zmiennych losowych \((X,X'), (Y,Y')\) o jednakowym rozkładzie i nie zakładamy ich niezależności. Wówczas \[Var\Big(\frac{g(X,Y)+g(X',Y')}{2}\Big)=\frac12Var(X,Y)\Big[1+corr(g(X,Y),g(X',Y'))\Big]=\frac12\sigma^2(1+\rho).\] Zatem, aby zmniejszyć wariancję, zmienne losowe muszą być ujemnie skorelowane.

W naszym przypadku, zamiast losować \(m\) zmiennych losowych \(X_i\sim U(0,1)\) i analogicznie dla \(Y\), możemy wylosować \(\frac m2\) zmiennych losowych \(X_i\sim U(0,1)\) i wyliczyć \(X'_i=1-X_i\) (analogicznie dla \(Y\)). Wtedy para \((X_i,X'_i)\) ma ten sam rozkład i jest ujemnie skorelowana.

Sprawdźmy, dla przykładu a), jak zachowuje się wartość całki i wariancja w zwykłej, i zmodyfikowanej metodzie MC.

# Porównanie zwyklej metody Monte Carlo z zmodyfikowana wersja - zmienne antytetyczne

f <- function(x,y) 1/(2*pi)*exp(-0.5*(x^2+y^2))  #funkcja gestosci dwuwymiarowygo rozkladu normalnego gdy \rho=0
m <- c(100, 1000, 10000, 100000) #liczba wylosowanych punktów

a=0
b=1
w1 = {}
w2 = {}
war1 ={}
war2 = {}
for (i in 1:4) {
  #tradycyjna MC
  x= a+(b-a)*runif(m[i])
  y= a+(b-a)*runif(m[i])
  w1[i] = (b-a)^2*mean(f(x,y))
  war1[i] = var(f(x,y))
  
  #MC ze zmiennymi antytetycznymi
  x1= runif(m[i]/2)
  y1= runif(m[i]/2)
  x2=1-x1
  y2=1-y1
  w2[i] = (b-a)^2*mean((f(x1,y1)+f(x2,y2))/2)
  war2[i] = var((f(x1,y1)+f(x2,y2))/2)
}
wartosc = 0.116516 # najlepsze przybliżenie wartosci calki
blad1 = abs(w1-wartosc)/wartosc #blad wzgledny MC
blad2 = abs(w2-wartosc)/wartosc #blad wzgledny MC Zm. antytetyczne
plot(c(2,3,4,5), blad1, type = 'o', col = 'red',ylim = c(0,0.025), 
     xlab = '10^n losowań', ylab = 'blad wzgledny przybliżenia', main = 'Porównanie bledów metody MC i MC zm. antet. na kwadracie')
lines(c(2,3,4,5),blad2,type = 'o', col= 'blue')

Czerwona krzywa - tradycyjna metoda MC; niebieska - metoda MC z zastosowaniem zmiennych antytetycznych.

podsumowanie <- data.frame("liczba losowań"=m,"wartosc MC"=w1, "wariancja MC"=war1, 
                           "blad wzgl MC"=blad1,"wartosc MC-z.antyt."=w2, "wariancja MC-z.antyt."=war2, 
                           "blad wzgl MC-z.antyt."=blad2)
podsumowanie
##   liczba.losowań wartosc.MC wariancja.MC blad.wzgl.MC wartosc.MC.z.antyt.
## 1          1e+02  0.1151677 0.0005391149 0.0115714654           0.1159442
## 2          1e+03  0.1156376 0.0005569758 0.0075387008           0.1163530
## 3          1e+04  0.1166705 0.0005547872 0.0013257016           0.1164738
## 4          1e+05  0.1164903 0.0005518396 0.0002201553           0.1165381
##   wariancja.MC.z.antyt. blad.wzgl.MC.z.antyt.
## 1          3.254720e-05          0.0049074248
## 2          2.534103e-05          0.0013991467
## 3          2.571408e-05          0.0003625357
## 4          2.583672e-05          0.0001900155

Jak można zauważyć, zmodyfikowana metoda Monte Carlo jest dokładniejsza i zmniejsza wariancję, w naszym przypadku, o rząd wielkości.

Źródło