Monte Carlo Integration

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

\[ E[g(x) = \int_{-\infty}^{\infty} g(x)f(x)dx\]

Hitunglah Estimasi Monte Carlo dari :

\[ \theta = \int_{0}^{1} e^{-x} dx \]

m<-10000

x<-runif(m)

theta.hat<-mean(exp(-x))

print(theta.hat)

## [1] 0.6309004

print(1-exp(-1))

## [1] 0.6321206

Untuk memberikan estimasi terhadap \[\int_{a}^{b} f(x)dx\] ada dua cara yang dapat dilakukan
Cara pertama melakukan transformasi fungsi sebagai akibat dari pengubahan limit integrasi dari a ke b menjadi 0 ke 1
Transformasi yang dilakukan adalah $$ \frac{x-a}{b-a} $$, maka\[ dy = \frac{dx}{b-a} \]

\[\int f(x) dx = \int_{0}^{1}f(y(b-a)(b-a)dy\]
Berikan estimator Monte Carlo dari \[\int_{4}^{2} e^{-x} dx\]
untuk \[f(x) = e^{-x} dx\], 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.1178093

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

## [1] 0.1170196

m<-10000

x<-runif(m,2,4)

theta.hat<-2*mean(exp(-x))

print(theta.hat)

## [1] 0.1169745

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

## [1] 0.1170196

Bahan Resume 1

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}{2}} sin x dx \]

monte_carlo <- function(f,a,b){
  m <- 10000
  y <- runif(m)
  theta.hat <- (b-a)* mean(f(y*(b-a)+a))
  return(theta.hat)
}

f <- function(x) 0.5*x^2
monte_carlo(f,2,6)

## [1] 34.66579

integrate(f,2,6)

## 34.66667 with absolute error < 3.8e-13

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

## [1] 0.5042391

integrate(f,0,pi/3)

## 0.5 with absolute error < 5.6e-15

Normal CDF

Generate bilangan acak Unifrom \[(0,1) u_1, u_2,..., U_m\]

Kemudian lakukan perhitungan

\[ \theta = \frac{}{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{\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.5398275 0.6430509 0.7366701 0.8157599 0.8780022 0.9236443 0.9548483
## 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.9747654 0.9866731 0.9933930
## Phi 0.9753892 0.9872365 0.9937903

Bahan Resume 2

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

Diketahui fungsi PDF

\[ f(x) = \int_0^\infty 5e^{-5t}dt \]

Hitunglah CDF dari distribusi exponensial tesebut untuk x=0 hingga x=5

\[ \theta = \int_0^e{-5t} dt \]
Lakukan Transformasi :

\[ y=\frac{t}{x}, t=xy \]

\[ \frac{dy}{dt}=\frac{1}{x} \]

\[ dt = x dy \]

Sehingga,

\[ \theta = \int_0^x e^{-5t} dt = \int_0^1 x e^{-5xy} dy \]

Maka, Estimasinya adalah 5 * 0

# 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.8355016 0.9793709 1.0063983 1.0132885 1.0164128 1.0186166
##  [8] 1.0204441 1.0220292 1.0234211 1.0246459 1.0257199 1.0266557 1.0274639
## [15] 1.0281540

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] = monte_carlo(f, 0, x[1])
  }
  else {
    p[i] = monte_carlo(f, x[i-1], x[i])
  }
  
}
cdf = cumsum(p)
cdf

##  [1] 0.0000000 0.8317964 0.9702998 0.9936701 0.9976167 0.9982731 0.9983835
##  [8] 0.9984020 0.9984051 0.9984056 0.9984057 0.9984057 0.9984057 0.9984057
## [15] 0.9984057

Resume KS Minggu 13

Monte Carlo Integration

Hitunglah Estimasi Monte Carlo dari :

Bahan Resume 1

Normal CDF

Bahan Resume 2