U statistici se često javljaju neki problemi koji nemaju analitička rešenja. Jedan od takvih problema je problem integracije. Za početak, posmatrajmo integral oblika \[ \int\limits_{\mathcal{X}}h(x)f(x)dx \] gdje je \(f\) neka funkcija gustine. Ovakav integral može se rešavati numeričkim metodama, sa kojima ste se već susreli u drugim kursevima.
U R-u postoje funkcije koje numeričkim metodama rešavaju integrale za funkciju koja se proslijedi, kao što su integrate() i area() (iz paketa MASS). Detaljnije o metodama koje kosriste i argumentima koji treba da se zadaju pročitajte u dokumentaciji R-a.

Primjeri

  1. Hoćemo da izračunamo vrijednost gama funckije funkcijom integrate (numeričkim metodama) i da je uporedimo sa pravim vrijednostima.
ch <- function(la) {
  integrate(function(x) {
    x ^ (la - 1) * exp(-x)
  }, 0, Inf)$val
} 
# Porede se logaritmovane vrijednosti 
plot(
  lgamma(seq(.01, 10, length.out = 100)),
  log(sapply(seq(.01, 10, length.out = 100), ch)),
  xlab = "log(integrate(f))",
  ylab = expression(log(Gamma(lambda))),
  pch = 19,
  cex = .6
)

  1. Posmatrajmo uzorak obima \(n=10\) iz Košijeve raspodjele sa parametrom \(\theta=350\). Hoćemo da posmatramo integral zajedničke gustine raspodjele po parametru \(\theta\) (ovakav neki problem može da se nađe u Bajesovom pristupu ocjenjivanju parametara o čemu ćete čuti na kursu Matematička statistika). \[ m(x)=\int\limits_{-\infty}^{\infty}\prod_{i=1}^{10}\frac{1}{\pi}\frac{1}{1+(x_i-\theta)^2}d\theta \]
cac <- rcauchy(10) + 350

lik <- function(theta) {
u <- dcauchy(cac, location = theta)
x <- prod(u)
return(x)
}

integrate(function(theta)
sapply(theta, lik), -Inf, Inf)
## 2.761064e-44 with absolute error < 5.5e-44
integrate(function(theta)
sapply(theta, lik), 200, 400)
## 8.144739e-11 with absolute error < 1.5e-10
# Ne dobijaju se dobre vrijednost(druga treba da je manja), a na to upucuje i izuzetno mala greska
# U ovom slucaju bolje ce raditi funkcija area
library(MASS)
cac<-rcauchy(10)
nin<-function(a){integrate(function(theta) sapply(theta,lik),-a,a)$val} #samo vrijednost integrala, bez dr podataka
nan<-function(a){area(function(theta) sapply(theta,lik),-a,a)} 
x<-seq(1,10^3,length.out=10^4)
y<-log(sapply(x,nin)) 
z<-log(sapply(x,nan))
plot(x,y,type="l",ylim=range(cbind(y,z)),lwd=2)
lines(x,z,lty=2,col="sienna",lwd=2)

MONTE KARLO METODE

“Hit or miss”

Ideja ove metode je da se simulira veliki broj tačaka (nezavisno) iz dvodimenzione uniformne raspodjele na nekom pravougaoniku na kojem se nalazi grafik funkcije koju želimo da integralimo. Kako je integral jednak površini ispod grafika funkcije, njegovu ocjenu ćemo dobiti tako što procenat tačaka koje se nađu ispod grafika pomnožimo sa površinom cijelog pravougaonika.
Mi ovdje zapravo posmatramo srednju vrijednost nezavisnih indikatora, pa će po zakonu velikih brojeva ona konvergirati kao očekivanju, s.s. Preciznije: Posmatramo integral \(\int\limits_{a}^{b}f(x)dx\). Pretpostavimo da se maksimum te funkcije postiže u nekoj tački \(c \in [a,b]\). Generišemo \((U_1^{i},U_2^{i}),\ i=1,2,\dots,n\), gdje je \(U_1 \in \mathcal{U}[a,b]\ \text{ i } \mathcal{U_2} \in \mathcal{U}[0,f(c)]\). Opisana ocjena će zapravo biti \(\frac{1}{n}\sum\limits_{i=1}^{n}I\{U_2\leq f(U_1)\}\) i po zakonu velikih brojeva ona konvergira ka \(E(I\{U_2 \leq f(U_1)\})=\frac{1}{(b-a)f(c)}\int\limits_{a}^{b}f(x)dx\), pa pomnoženo sa površinom pravougaonika \([a,b]\times[0,f(c)]\) je baš jednako integralu koji ocjenjujemo. \[ E(I\{U_2 \leq f(U_1)\})=E(EI\{U_2 \leq f(U_1)\}|U_1)=\frac{1}{b-a}\int\limits_{a}^{b}P\{U_2 \leq f(U_1)|U_1=u_1\}du_1=\frac{1}{f(c)(b-a)}\int\limits_{a}^{b}f(x)dx \]

Napomena: Kako po CGT ova suma indikatora ima normalnu graničnu raspodjelu, ako želimo da unaprijed zadamo neku grešku \(\varepsilon\), pomoću kvantila normalne raspodjele možemo da odredimo najmanji broj \(n\) takav da je sa vjerovatnoćom \(1-\alpha\) greška integracije manja od \(\varepsilon\). Pokažite da ako ograničimo disperziju indikatora sa \(\frac{1}{4}\) i posmatramo simetričan interval povjerenja za ocjenu integrala, dobija se da je to \(n=[(f(c)*(b-a)*z_{1-\alpha/2})^2/4\varepsilon^2]+1\). Ja ću u primjerima najčešće koristiti proizvoljno izabrano \(n\).

Primjeri

  1. Izračunati integral \(\int\limits_{0}^{2}e^{-\frac{x^2}{2}}dx\).
curve(exp(-(x ^ 2 / 2)), 0, 2)

# max je 1 pa tražimo slučajne tačke na pravougaoniku [0,2]x[0,1]
# neka je d=f(c), tj max 
mk.integral <- function(n, a, b, d, f) {
  u1 <- runif(n, a, b)
  u2 <- runif(n, 0, d)
  mean(u2 <= f(u1)) * (b - a) * d
  
}

mk.integral(10000, 0, 2, 1, function(x)
exp(-(x ^ 2 / 2)))
## [1] 1.2026
integrate(function(x)
exp(-(x ^ 2 / 2)), 0, 2)
## 1.196288 with absolute error < 1.3e-14
# ili
sqrt(2 * pi) * (pnorm(2) - pnorm(0))
## [1] 1.196288
  1. Aproksimacija broja \(\pi\).
curve(sqrt(1-x^2),0,1)

# Iskoristimo funckiju iz prethodnog primjera
4 * mk.integral(100000, 0, 1, 1, function(x)
  sqrt(1 - x ^ 2))
## [1] 3.1372

Monte Karlo integracija

Vratimo se na problem računanja očekivanja neke slučajne veličine \[ E(h(X))=\int\limits_{\mathcal{X}}h(x)f(x)dx \] Ako generišemo uzorak iz te raspodjele koju ima \(X\) i uzmemo uzoračku sredini, po zakonu velikih brojeva, ona će konvergirati skoro sigurno ka integralu \(E(h(X))\).

  1. Za isti primjer \(\int\limits_{0}^{2}e^{-\frac{x^2}{2}}dx\) možemo uzeti da je gustina \(f(x)=1/2,\ \text{na } [0,2]\), a \(h(x)=2e^{-\frac{x^2}{2}}\), tj. treba da generišemo uzorak iz uniformne \(\mathcal{U}[0,2]\) raspodjele.
u <- runif(10000, 0, 2)
x <- 2 * exp(-(u ^ 2 / 2))
mean(x)
## [1] 1.18107
  1. Računanje funkcije standardne normalne raspodjele u nekoj tački t.

\[ \Phi(t)=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{t}e^{-\frac{x^2}{2}}dx=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{\infty}I\{x \leq t\}e^{-\frac{x^2}{2}}dx \] Generisaćemo brojeve iz normalne raspodjele i tražiti očekivanje indikatora iz izraza gore.

u<-rnorm(10^7)
g<-function(t){
  u<=t
}
mean(g(0))
## [1] 0.5000795
  1. Izračunati \(P\{X>20\}\) ako je \(X \in \mathcal{N}(0,1)\). Kako je ovo vrijednost veoma mala, ako iskoristimo prethodnu metodu dobićemo uvijek vrijednost 0.
1-mean(g(20))
## [1] 0

Kako je normalna raspodjela simetrična, \(P\{X>20\}=P\{X<-20\}\).Zato ćemo posmatrati intgral oblika \(\int\limits_{-\infty}^{t}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}dx,\ t<0\) koji se smjenom \(x=\frac{1}{u}\) svodi na sledeće: \[ \int\limits_{\frac{1}{t}}^{0}\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2u^2}}\frac{1}{u^2}du \] Nama u zadatku treba specijalno za \(t=-20\). Sada se ovo svodi na očekivanje \(g(U)\) gdje je je \(U \in \mathcal{U}[-\frac{1}{20},0]\)

u<-runif(10^5,-1/20,0)
g<-function(t){
  1/sqrt(2*pi)*exp(-1/(2*(t^2)))*1/(t^2)
}
mean(1/20*g(u))
## [1] 2.931518e-89

Importance sampling

“Uzorkovanje po značajnosti” je tehnika uz pomoću koje se smanjuje disperzija greške nastale prilikom ocjene srednje vrijednosti Monte Karlo metodom. Klasična metoda za ocjenjivanje srednje vrijednosti E(g(X)) za neku realnu funkciju \(g\) i slučajnu veličinu \(X\), iz uzorka \((X_1,...,X_n)\) bila bi: \(\frac{1}{n}(g(X_1)+...+g(X_n))\). Kako su neke vrijednosti slučajne veličine \(X\) značajnije za parametar koji se ocjenjuje, u ovom slučaju \(E(g(x))\), ideja je da se pronđe raspodjela koja ce takvim vrijednostima da dodjeli najveće vjerovatnoće i da se preko takve raspodjele ocijeni traženi parametar. Ispostavlja se da će se time sto ce se “značajnije’ vrijednosti pojavljivati češće u uzorku, smanjiti disperzija ocjene.
Ova tehnika je zasnovana na sledećim rezultatima:
T1. Neka je \(X\) diskretna slučajna veličina i \(g\) neka realna funkcija. Neka je Y neka druga slučajna velicina takva da važi \(P\{Y=x\}>0\), kad god je \(P\{X=x\}>0\). Označimo \(p_X(x)=P\{X=x\}\). Tada: \[ E(g(X)) = E[\frac{p_X(Y)}{p_Y(Y)}g(Y)] \] T2. Neka je \(X\) apsolutno neprekidna slučajna velicina i \(g\) neka realna funkcija. Neka je \(Y\) neka druga slučajna veličina takva da važi \(f_Y(x)>0\), kad god je \(f_X(x)>0\). Tada: \[ E(g(X)) = E[\frac{f_X(Y)}{ f_Y(Y)}g(Y)] \] Primjer: Neka je \(X \in \mathcal{N}(0,1)\) i \(g(x)=10e^{-5(x-3)^4}\). Ocijeniti srednju vrijednost \(E(g(X))\) klasinom Monte Karlo metodom i metodom Importance sampling.

# Resenje: Funckija g(x) dostize svoj maksimum u tacki x=3 i potom naglo opada 
# ka nuli za sve xeve manje i vece od 3. S druge strane, gustina normalne
# raspodjele u tacki x=3 je veoma blizu nuli. Dakle, ako budemo ocjenjivali
# ocekivanje od g(X) standardnom Monte Karlo metodom, imacemo veliki broj
# vrijednosti iz uzorka za koje ce g(x) biti skoro 0, a mali broj vrijednosti za
# koje ce ona biti razlicita od nule sto za rezultat ima veliku disperziju
# ocjene.

n<-10000
g<-function(x){
  exp(-5*(x-3)^4)
}

x <- rnorm(n)      # uzorak iz N(0,1)
mk<-mean(g(x))     # Monte Karlo ocjena srednje vrijednosti
d1<-1/n*var(g(x))  # disperzija ocjene mk
sd1<-sqrt(d1)      # standardno odstupanje ocjene mk

# Kako se maksimum fje g(x) postize u tacki x=3, za y uzimamo normalnu
# raspodjelu sa ocekivanjem 3, sto je ujedno i moda te raspodjele (vrijednost sa
# najvecom vjerovatnocom), pa odgovara ideji za importance sampling

y<-rnorm(n, mean = 3)
is<-mean(dnorm(y,0,1)/dnorm(y,3,1)*g(y))
d2<-1/n*var(dnorm(y,0,1)/dnorm(y,3,1)*g(y))
sd2<-sqrt(d2)

# Poredimo standardna odstupanja ove dvije ocjene:

sd1-sd2
## [1] 0.0006488373