У статистици се често јављају неки проблеми који немају аналитичка решења. Један од таквих проблема је проблем интеграције. Посматрајмо интеграл облика \[\int_X h(x)f(x) dx ,\] где је \(f(x)\) нека функција густине. Овакав интеграл се може решавати нумеричким методама.

У R-у постоје функције које нумеричким методама решају интеграле за функцију која им се проследи, као што су integrate() и area() (из пакета MASS).

Примери

ch <- function(la) {
  integrate(function(x){x^(la-1)*exp(-x)}, 0, Inf)$val # vraca ocenjene vrednosti integrala
} 
# Poredimo graficki
plot(
  gamma(seq(.01, 10, length.out = 100)),
  sapply(seq(.01, 10, length.out = 100), ch),
  xlab = "integrate(f)",
  ylab = expression(Gamma(lambda)),
  pch = 19,
  cex = .6
)

Монте Карло методе

Oве нумеричке методе код којих се употребљавају низови случајних бројева за извршење симулације. Првобитно су познати као ”статистичка упрошћавања”, али назив ”Монте Карло”, популаризован од стране првих истраживача у овој области је проистекао из назива чувеног казина у Монаку.

Методама Монте Карло се могу решити и неки задаци у којима се класичне методе нумеричке математике не могу применити.

Оно што је такође значајно за алгоритме Монте Карло јесте да су обично једноставни и лаки за програмирање.

Интеграција Монте Карло

Интеграција Монте Карло је Монте Карло метода којом се (приближно) нумерички израчунава дати интеграл. Примењује се када је интеграл сложен и аналитички тежак или немогућ за израчунавање.

За приближно израчунавање интеграла \[I = \int_{a}^{b}g(x) dx\] користе се две Монте Карло методе:

Монте Карло метода погодака и промашаја (The Hit and Miss Monte Carlo Method)

Нека је функција \(g(x)\) ограничена \(0 \leq g(x) \leq c\), \(a \leq x \leq b\) и означимо са \(D\) област \[D=\{(x,y):a \leq x \leq b, 0 \leq y \leq c\}.\]

Нека је \((X,Y)\) случајни вектор са униформном расподелом на правоугаонику \(D\) са густином расподеле

\[ f_{X,Y}(x,y) = \frac{1}{c(b-a)}, \quad (x,y) \in D\]

Означимо са \(S\) површ испод криве \(g(x)\), \(S=\{(x,y):y \leq g(x)\}\). Вероватноћа да се случајни вектор нађе испод криве \(g(x)\) је тада једнака \[p=\frac{\int_{a}^{b} g(x) dx}{c(b-a)} = \frac{I}{c(b-a)}\]

Ако претпоставимо да је генерисано \(N\) независних случајних вектора (\(X_1,Y_1), (X_2,Y_2),...,(X_N,Y_N)\), на основу закона великих бројева, вероватноћа \(p\) се може оценити са: \[\hat{p} = \frac{N_H}{N},\] где је:

  • \(N_H\) број погодака тј.број случајева када је \(Y_i \leq g(X_i),\quad i = 1,2,…,N\)

  • \(N-N_H\) број промашаја тј. број случајева када је \(Y_i > g(X_i), \quad i = 1,2,…,N\).

Из претходног закључујемо да се интеграл \(I\) може оценити као: \[I \approx \theta = c(b-a)\frac{N_H}{N},\]

Сваки од \(N\) покушаја има Бернулијеву расподелу са вероватноћом \(p\), па случајна величина \(N_H\) има Биномну расподелу са параметрима \((N,p)\).

Математичко очекивање и дисперзија оцене \(\theta\) су: \[E(\theta) = \frac{c(b-a)}{N}E(N_H)= c(b-a)p = I\] \[D(\theta) = \frac{c^2(b-a)^2}{N^2}D(N_H) = c^2(b-a)^2 \frac{p(1-p)}{N}\] Треба приметити да је оцена \(\theta\) непристрасна, а пошто \(D(\theta)\rightarrow0\), када \(N \rightarrow \infty\), оцена \(\theta\) је и постојана.

Алгоритам за Монте Карло методу погодака и промашаја

Приметимо да је \(N_H = \sum_{i=1}^n I\{ Y_i \le g(X_i)\}\), па применом ЦГТ ова сума индикатора има граничну нормалну расподелу. Ако желимо унапред да задамо неку грешку \(\epsilon\), помоћу квантила нормалне расподеле можемо да одредимо најмањи број \(n\) такав да је са вероватноћом \(1-\alpha\) грешка интеграције мања од \(\epsilon\).

  • Ако претпоставимо да је дисперзија индикатора мања од \(\frac{1}{4}\) и да је интервал поверења за оцену интеграла симетричан, добија се да је \(N=\left[\frac{(g(c)(b-a)z)^2}{4\epsilon^2}\right]+1\).

  • Генерисати низ \((U_j)_{j=1}^{2N}\) од \(2N\) псеудослучајних бројева.

  • Поређати бројеве у \(N\) парова \((U_1,U_1^{'}), (U_2,U_2^{'}),...,(U_N, U_N^{'})\) на било који начин тако да се сваки број \(U_i\) појави тачно једанпут.

  • Израчунати \(X_i=a+(b-a)U_i\), \(g(X_i)\) и \(Y_i = cU_i^{'}\) за \(i=1,2,...N.\)

  • Пребројати број погодака \(N_H\) за које важи \(Y_i \leq g(X_i)\).

  • Оценити интеграл \(I\) са \[\theta = c(b-a)\frac{N_H}{N}\]

Монте Карло метода погодака и промашаја - пример

Користећи Монте Карло методу погодака и промашаја израчунати интеграл \[I = \int_{0}^{2} e^{-x^2} dx.\] За грешку апроксимације узети \(\epsilon = 0,001\) и ниво поверења \(\alpha = 0,95\).

integral_e1 <- function(){
c<-1 # jer je podintegralna funkcija manja od 1
a<-0
b<-2 # granice integrala
epsilon <- 0.001 # greska aproksimacije
alfa <- 0.95 # nivo poverenja
z <- qnorm((1+alfa)/2)
N <- trunc(((c*(b-a)*z)^2)/(4*epsilon^2))+1 # obim uzorka
U <- runif(N) #generisemo niz od N clanova
V <- runif(N) #generisemo niz od N clanova
X <- a+(b-a)*U
Y <- c*V
Nh <- sum(Y<=exp(-X^2))
teta1 <- c*(b-a)*Nh/N
return (teta1)
}
integral_e1()
## [1] 0.8820685

Монте Карло метода узорачке средине (The Sample Mean Monte Carlo Method)

Број потребних експеримената \(N\) се може одредити тако да је: \[P\{|\theta - I|\leq \epsilon\} = \alpha,\] где је:

  • \(\epsilon\) – грешка апроксимације,

  • \(\alpha\) – ниво поверења, (обично се узима \(\alpha = 95\%\) или \(\alpha = 99\%\)).

За довољно велико \(N\), применoм централнe граничнe теоремe имамо: \[\theta^{*} = \frac{\theta-I}{\sigma_{\theta}} \sim \mathbf{N}(0,1).\] Тада се добија \(P\{|\theta^{*}| \leq z_{\alpha}\}=\alpha\) где \(z_{\alpha}\) налазимо тако да је \(F(z_{\alpha}) = \frac{1+ \alpha }{2}\), где је \(F\) функција нормалне расподеле.

Користећи апроксимацију \(p(1-p)\approx \frac{1}{4}\), добијамо: \[N \geq \frac{c^2(b-a)^2z_{\alpha}^2}{4\epsilon^2}.\]

Интеграл \(I = \int_{a}^{b} g(x) dx\) се може представити као очекивана вредност неке случајне величине. Прво напишимо интеграл на следећи начин: \[I = \int_{a}^{b} \frac{g(x)}{f_X(x)} f_X(x) dx.\] Затим претпоставимо да је \(f_X(x)\) произвољна густина расподеле таква да је \(f_X(x) >0\) када \(g(x)\neq 0\). Тада је: \[I = E\left(\frac{g(X)}{f_X(X)}\right),\] где је \(X\) случајна величина са густином расподеле \(f_{X}(x)\).

Узмимо, једноставности ради, да \[ f_{X}(x) = \frac{1}{b-a}, \quad x \in (a,b).\]

Тада имамо да је \(E(g(X)) = \frac{I}{b-a}\), односно \[I = (b-a)E(g(X)).\] Oцена за интеграл \(I\) је узорачка средина \[\theta = (b-a) \cdot \frac{1}{N} \sum_{i=1}^{N} g(X_i).\]

Математичко очекивање и дисперзија оцене \(\theta\) су: \[E(\theta) = I\] \[D(\theta) = ... = \frac{1}{N} \left( (b-a) \int_{a}^{b} g^2(x) dx - I^2 \right)\] Треба приметити да је оцена \(\theta\) непристрасна, а пошто \(D(\theta)\rightarrow0\), када \(N \rightarrow \infty\), оцена \(\theta\) је и постојана.

Алгоритам за Монте Карло методу узорачке средине

  • Генерисати низ \((U_j)_{j=1}^N\) од \(N\) псеудослучајних бројева.

  • Израчунати \(X_i = a+(b-a)U_i\), \(i=1,2,...,N\).

  • Израчунати \(g(X_i)\), \(i=1,2,...,N\).

  • Оценити интеграл узорачком средином \[\theta = (b-a) \frac{1}{N} \sum_{i=1}^{N} g(X_i).\]

Монте Карло метода узорачке средине - пример

Користећи Монте Карло методу узорачке средине израчунати интеграл \[I = \int_{0}^{2} e^{-x^2} dx.\]

integral_e2 <- function(){
a<-0
b<-2 # granice integrala
N <- 10000 #obim uzorka
U <- runif(N) #generisemo N slucajnih brojeva
X <- a+(b-a)*U
teta2 <- (b-a)*1/N*sum(exp(-X^2))
return(teta2)
}
integral_e2()
## [1] 0.8740199

Апроксимација броја \(\pi\)

Ако насумице изаберемо тачку у равни унутар јединичног квадрата \([0,1]X[0,1]\), вероватноћа да ће се она наћи унутар јединичне кружнице је однос површине четвртине круга и укупне површине квадрата. \[p = \frac{P_{deoKruga}}{P_{kvadrata}} = \frac{\frac{1^2\pi}{4}}{1^2} = \frac{\pi}{4}\]

Закон великих бројева нам говори да ће однос броја тачака које су завршиле унутар јединичне кружнице и укупног броја тачака бити приближно једнака траженој вероватноћи.