Monte Carlo Integration

Integral Monte Carlo adalah integral menggunakan metode statistika berdasarkan sampel acak.
Misalkan \(X\) adalah peubah acak dengan fungsi kepadatan \(f(x)\)
Ekspektasi dari peubah acak \(Y = g(X)\) adalah:

\(E[g(X)] = \int_{-\propto}^{\propto} g(x)f(x)dx\)
Pandang permasalahan estimasi integral \(\theta = \int_{0}^{1} g(x)dx\)
Jika \(X_{1},X_{2},\cdots,X_{m}\) adalah sampel acak Uniform(0,1), maka

\(\widehat{\theta} = \overline{g_{m}(X)} = \frac{1}{m}\sum_{i=1}^{m}g(X_{i})\)

konvergen ke \(E[h(X)]\)
Estimator Monte Carlo dari \(\int_{0}^{1}g(x)dx\) adalah \(\frac{1}{m}\sum_{i=1}^{m}g(X_{i})\)
Hitunglah estimasi Monte Carlo dari

\(\theta = \int_{0}^{1}e^{-x}dx\)

dan bandingkan hasilnya dengan nilai eksak.

m<-10000
y<-runif(m)
theta.hat<-2*mean(exp(-2*y-2))
print(theta.hat)

## [1] 0.1172285

print(exp(-2)-exp(-4))

## [1] 0.1170196

m = 2
x = runif(m)
theta.hat = mean(exp(-x))
print(theta.hat)

## [1] 0.5665544

print(1-exp(-1))

## [1] 0.6321206

Untuk memberikan estimasi terhadap \(\int_{a}^{b}f(x)dx\), ada dua cara yang dapat dilakukam.
Cara pertama, melakukan transformasi fungsi sebagai akibat dari pengubahan limit integrasi dari a ke b menjadi 0 dan 1.
Transformasi yang dilakukan adalah \(y\frac{x-a}{b-a}\). Maka \(dy=\frac{dx}{b-a}\) \(\int_{a}^{b}f(x)dx=\int_{0}^{1}f(y(b-a)+1)(b-a)dy\)
Berikan estimator Monte Carlo dari \(\int_{2}^{4}e^{-x}dx\)
Untuk \(f(x)=e^{-x}\), \(a=2\), dan \(b=4\), maka \(b-a=2\).
Sehingga \(f(y(b-a)+a)(b-a)=2f(2y+2)=2e^{-2y-2}\) \(\int_{2}^{4}e^{-x}dx=\int_{0}^{1}2e^{-2y-2}dy\)

m<-10000
y<-runif(m)
theta.hat<-2*mean(exp(-2*y-2))
print(theta.hat)

## [1] 0.1163675

print(exp(-2)-exp(-4))

## [1] 0.1170196

Cara kedua adalah dengan mengganti kepadatan Uniform (0,1) menjadi kepadatan Unifrom (a,b)
Namun dengan mengubah kepadatan, artinya fungsi awal dikalikan dengan \(b-a:\widehat{\theta}=(b-a)\overline{g(X)}\)

m<-10000
x<-runif(m,2,4)
theta.hat<-2*mean(exp(-x))
print(theta.hat)

## [1] 0.1168497

print(exp(-2)-exp(-4))

## [1] 0.1170196

Bahan resume: Buatlah fungsi dasar integrasi Monte Carlo dengan argumen fungsi f, batas bawah a, dan batas atas b.

Dari fungsi Monte Carlo yang sudah anda buat, berikan estimasi dari integral berikut: - \(\int_{2}^{6}\frac{dx}{x^{2}}\) - \(\int_{0}^{\frac{\pi}{3}}sin x dx\) fungsi Monte Carlo

MonteCarlo <- function(f,a,b){
  m <- 200000
  x <- runif(m,a,b)
  theta <- mean((b-a)*f(x))
  print(theta)
}
f <- function(x) 1/(x^2)
MonteCarlo(f,0,1)

## [1] 152989.9

f<-function(x) sin(x)
MonteCarlo(f,0,pi/3)

## [1] 0.500863

Gunakan pendekatan Monte Carlo untuk memberikan estimasi fungsi kepadatan normal baku berikut:

$F(x)=\int_{-\propto}^{x}\frac{1}{\sqrt{2\pi}}e^{\frac{-t^{2}}{2}}$

Bagi kasus untuk \(x \ge 0\) dan \(x \lt 0\). Karena distribusi normal baku, maka untuk \(x \ge 0\) nilai cdf-nya adalah 0.5
Cukup memberikan estimasi untuk \(\theta = \int_{0}^{x}e^{\frac{-t^{2}}{2}}dt\) untuk \(x \ge 0\).
Lakukan transformasi fungsi dengan mengikuti pengubahan batas dari 0 ke x menjadi 0 ke 1.
Dapat digunakan transformasi \(y = \frac{t}{x}\). Sehingga diperolah \(dt = x dy\).

\(\theta = \int_{0}^{x}e^{\frac{-t^{2}}{2}}dt = \int_{0}^{1}xe^{\frac{-(xy)^{2}}{2}}dy\)
Generate bilangan aca Uniform(0,1)\(u_1,u_2,\dots,u_m\)
Kemudian lakukan perhitungan

\(\widehat{\theta}=\overline{g_m(u)}=\frac{1}{m}\sum_{i=1}^{m}xe^{\frac{-(u_ix)^{2}}{2}}\)
Maka estimasi dari Monte Carlo adalah \(F(x) = 0.5 + \frac{\widehat{\theta}}{\sqrt{2\pi}}\)

x<-seq(.1,2.5,length=10)
m<-10000
u<-runif(m)
cdf<-numeric(length(x))
for(i in 1:length(x)){
  g<-x[i]*exp(-(u*x[i])^2/2)
  cdf[i]<-0.5+mean(g)/sqrt(2*pi)
}

Sudah didapatkan nilai estimasi Monte Carlo dari cdf.
Dalam R terdapat juga fungsi dasar yang menghitung cdf pada titik tertentu, yaitu pnorm.
Bandingkan hasil estimasi Monte Carlo dengan pnorm

Phi<-pnorm(x)
print(rbind(x,cdf,Phi))

##          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
## x   0.1000000 0.3666667 0.6333333 0.9000000 1.1666667 1.4333333 1.7000000
## cdf 0.5398269 0.6430230 0.7365294 0.8153689 0.8771869 0.9222135 0.9526144
## Phi 0.5398278 0.6430662 0.7367420 0.8159399 0.8783275 0.9241187 0.9554345
##          [,8]      [,9]     [,10]
## x   1.9666667 2.2333333 2.5000000
## cdf 0.9715631 0.9823731 0.9879095
## Phi 0.9753892 0.9872365 0.9937903

Bahan Resume:

Cobakan untuk fungsi yang berbeda dan nilai x yang berbeda.
Bagaimanakah jika tidak dilakukan transformasi fungsi? Yaitu ketika batasnya tetap dari 0 ke x. Ubah kodenya dalam R.

# CDF distribusi eksponensial dengan monte carlo
x<-seq(0 , 5,length=15)
m<-10000
u<-runif(m)
cdf<-numeric(length(x))
for(i in 1:length(x)){
  g = x[i] * exp(-5*u*x[i])
  cdf[i]<- mean(g) * 5
}

cdf

##  [1] 0.0000000 0.8322470 0.9696410 0.9900936 0.9914176 0.9901360 0.9889389
##  [8] 0.9881699 0.9877952 0.9877374 0.9879284 0.9883126 0.9888447 0.9894876
## [15] 0.9902103

pexp(x, 5)

##  [1] 0.0000000 0.8323228 0.9718843 0.9952856 0.9992095 0.9998675 0.9999778
##  [8] 0.9999963 0.9999994 0.9999999 1.0000000 1.0000000 1.0000000 1.0000000
## [15] 1.0000000

# Tanpa Transformasi
x<-seq(0 , 5,length=15)
f = function(x) 5*exp(-5*x)
p <-numeric(length(x))
for(i in 1:length(x)){
  if (i == 1){
    p[i] = MonteCarlo(f, 0, x[1])
  }
  else {
    p[i] = MonteCarlo(f, x[i-1], x[i])
  }
  
}

## [1] 0
## [1] 0.8338708
## [1] 0.139395
## [1] 0.02337229
## [1] 0.003916239
## [1] 0.0006585125
## [1] 0.0001104219
## [1] 1.848415e-05
## [1] 3.100751e-06
## [1] 5.200837e-07
## [1] 8.719071e-08
## [1] 1.46393e-08
## [1] 2.452926e-09
## [1] 4.110716e-10
## [1] 6.881718e-11

cdf = cumsum(p)
cdf

##  [1] 0.0000000 0.8338708 0.9732658 0.9966381 1.0005543 1.0012128 1.0013232
##  [8] 1.0013417 1.0013448 1.0013453 1.0013454 1.0013454 1.0013454 1.0013454
## [15] 1.0013454

Monte Carlo Integration

Komputasi Statistik - Sains Data ITERA

TERIMA KASIH