Pengantar Metode Monte Carlo

Pengantar

Dua kelas utama masalah numerik yang muncul dalam prosedur analisis data adalah masalah optimisasi dan integrasi. Tidak selalu memungkinkan untuk menghitung secara analitik penduga yang terkait dengan suatu model, dan kita sering beralih ke solusi numerik. Salah satu cara untuk mengatasi masalah tersebut adalah dengan menggunakan simulasi. Estimasi Monte Carlo merujuk pada mensimulasikan pengambilan sampel hipotetis dari suatu distribusi probabilitas untuk menghitung besaran-besaran penting dari distribusi tersebut.

Gagasan dasar Monte Carlo adalah menuliskan integral sebagai nilai harapan terhadap suatu distribusi probabilitas, lalu mendekatinya menggunakan penduga metode momen ($E[g(X)] \approx \overline{g(X)} = \dfrac{1}{n}\sum g(X_{i})$).

Jika kita memiliki fungsi kontinu $g(\theta)$ dan kita ingin mengintegralkannya pada interval (a,b), kita dapat menuliskan ulang integral sebagai nilai harapan dari distribusi uniform $U \sim U[a,b]$, yaitu:

Menggunakan penduga metode momen, pendekatan integral kita adalah:

Dengan

adalah nilai-nilai hasil simulasi dari distribusi uniform.

Contoh 1: Pendekatan integral eksponensial

1) Diberikan fungsi $f(x)= e^{x}$, integral pada interval [3,5] adalah:

Pendekatan Monte Carlo untuk integral tersebut adalah:

#       Declaring the desired function
f = function(x){return(exp(x))}
#       Declaring the absolute error function
error = function(x,y){return(abs(x-y))}
#       The actual integral answer
ans = exp(5)-exp(3)

set.seed(6971)
#       number of iterations
n = 10^2
#       simulated uniform data
x= runif(n,3,5)
#       MonteCarlo approximation
MCa= (5-3)*mean(f(x))
#       Approximation error
e = error(ans,MCa)

rest =  data.frame(n = n,MCapprox = MCa,error = e)
set.seed(6971)
for(k in 3:6){
  n = 10^k
  x = runif(n,3,5)
  mca = (5-3)*mean(f(x))
  rest = rbind(rest,c(n,mca,error(ans,mca) ) )
}

kable(rest,digits = 5,align = 'c',caption = "Integral Monte Carlo approximation results",
      col.names =c("Number of simulations","Monte Carlo approximation","Error approximation"))Didukung oleh

Pendekatan Monte Carlo yang digeneralisasi

Secara umum, pendekatan integral untuk suatu distribusi f adalah:

Algoritme untuk membangun $\widehat{I}$ dapat dijelaskan dengan langkah-langkah berikut:

1) Hasilkan dari suatu distribusi f

2) Hitung:

3) Peroleh rataan sampel: $$\overline{I} = \dfrac{1}{n}\sum_{k=1}^{n}\dfrac{g(\theta_k)}{f(\theta_k)}$$

Pada bagian berikut, fungsi pendekatan Monte Carlo sederhana disajikan untuk menunjukkan cara kerja algoritme, dengan a dan b sebagai parameter densitas uniform, n jumlah simulasi yang diinginkan, dan f adalah fungsi yang ingin kita integralkan.

# The simple Monte Carlo function
MCaf = function(n,a,b,f){
x = runif(n,a,b)
MCa = (b-a)*mean(f(x))
return(MCa)
}Didukung oleh

Metode Monte Carlo dalam Analisis Data Bayesian

Gagasan utama analisis data Bayesian adalah menyesuaikan model ( seperti model regresi atau deret waktu ) menggunakan pendekatan inferensi Bayesian. Kita mengasumsikan bahwa parameter yang diminati memiliki distribusi teoretis; distribusi ini (posterior) diperbarui menggunakan distribusi data teramati (likelihood) dan informasi sebelumnya atau eksternal tentang parameter kita (prior distribution) dengan menggunakan teorema Bayes.

$$p(\theta / X) \text{ } \alpha\text{ } p(X/\theta)p(\theta)$$ Di mana:

$p(\theta/X)$ adalah distribusi posterior dari parameter

$p(X/\theta)$ adalah distribusi penarikan sampel dari data teramati ( likelihood ).

$p(\theta)$ adalah distribusi prior dari parameter.

Masalah utama dalam pendekatan Bayesian adalah mengestimasi distribusi posterior. Metode Markov Chain Monte Carlo ( mcmc ) menghasilkan sampel dari distribusi posterior dan mendekatkan nilai harapan, probabilitas, atau kuantil menggunakan metode Monte Carlo.

Pada dua bagian berikut, kami memberikan dua contoh untuk mendekati probabilitas dan kuantil dari suatu distribusi teoretis. Prosedur yang disajikan di atas merupakan metodologi lazim yang digunakan dalam pendekatan Bayesian.

Contoh 2: Pendekatan probabilitas dari distribusi gamma

Misalkan kita ingin menghitung probabilitas bahwa peubah acak $\theta$ berada antara nol dan 5 $P(0 < \theta < 5)$, dengan $\theta$ berdistribusi gamma dengan parameter a = 2 dan b = 1/3 ($\theta \sim Gamma(a = 2,b = 1/3)$), maka probabilitasnya adalah:

Dengan $I_{[0,5]}(\theta) = 1$ jika $\theta$ berada pada interval [0,5]. Gagasan pendekatan Monte Carlo adalah menghitung jumlah observasi yang berada pada interval [0,5], kemudian membaginya dengan total data hasil simulasi.

set.seed(6972)
#       number of iterations
n = 10^2
#       simulated uniform data
x= rgamma(n,shape = 2,1/3)
#       MonteCarlo approximation
MCa= mean(x <= 5)
#       Approximation error
e = error(pgamma(5,2,1/3),MCa)

rest =  data.frame(n = n,MCapprox = MCa,error = e)

for(k in 3:6){
  n = 10^k
  x= rgamma(n,shape = 2,1/3)
  mca= mean(x <= 5)
 rest = rbind(rest,c(n,mca,error(pgamma(5,2,1/3),mca) ) )
}

kable(rest,digits = 5,align = 'c',caption = "Probability Monte Carlo approximation results",
      col.names =c("Number of simulations","Monte Carlo approximation","Error approximation"))Didukung oleh

Contoh 3: Pendekatan kuantil dari distribusi normal

Misalkan kita ingin menghitung kuantil 0,95 dari peubah acak $\theta$ yang berdistribusi normal dengan parameter $\mu = 20$ dan $\sigma = 3$ ($\theta \sim normal(\mu = 20,\sigma^{2} = 9)$), maka kuantil 0,95 adalah:

Gagasan utamanya adalah menemukan nilai sampel terbesar yang memberikan probabilitas sama dengan atau kurang dari 0,95; pendekatan kuantil Monte Carlo memperkirakannya menggunakan fungsi quantile() atas data hasil simulasi.

set.seed(6973)
#       number of iterations
n = 10^2
#       simulated uniform data
  x= rnorm(n,20,3)
#       MonteCarlo approximation
MCa= quantile(x,0.95)
#       Approximation error
e = error(qnorm(0.95,20,3),MCa)

rest =  data.frame(n = n,MCapprox = MCa,error = e)

for(k in 3:6){
  n = 10^k
  x= rnorm(n,20,3)
  mca= quantile(x,0.95)
 rest = rbind(rest,c(n,mca,error(qnorm(0.95,20,3),mca) ) )
}

kable(rest,digits = 5,align = 'c',caption = "Quantile Monte Carlo approximation results",
      col.names =c("Number of simulations","Monte Carlo approximation","Error approximation"),row.names = FALSE)Didukung oleh

Diskusi dan Kesimpulan

Metode pendekatan Monte Carlo menawarkan alat alternatif untuk pendekatan integral dan merupakan alat penting dalam pendekatan inferensi Bayesian, terutama ketika kita bekerja dengan model yang canggih dan kompleks. Seperti terlihat pada ketiga contoh kita, metode Monte Carlo memberikan pendekatan yang sangat baik, namun membutuhkan jumlah simulasi yang sangat besar untuk mendapatkan galat pendekatan yang mendekati nol.

Referensi

1) Introducing Monte Carlo methods with R, Springer 2004, Christian P. Robert and George Casella.

2) Handbook of Markov Chain Monte Carlo, Chapman and Hall, Steve Brooks, Andrew Gelman, Galin L. Jones, and Xiao-Li Meng.

3) Introduction to mathematical Statistics, Pearson, Robert V. Hogg, Joseph W. Mckean, and Allen T. Craig.

4) Statistical Inference An Integrated Approach, Chapman and Hall, Helio S. Migon, Dani Gamerman, Francisco Louzada.