Aproksimasi Monte Carlo ~ taudata Academy

Apakah metode Monte Carlo itu? Bagaimana cara kerjanya? Pada situasi seperti apa saya memerlukan metode Monte Carlo? Untuk menjawab pertanyaan - pertanyaan tersebut, perhatikan contoh berikut. Pada contoh model Poisson (link), diperoleh distribusi posterior dari parameter $\theta_A$ dan $\theta_B$ adalah:

$\theta_A|n_A=40, \sum{Y_{i,A}}=36\sim dgamma(\theta_A, 38,43)$

$\theta_B|n_B=25, \sum{Y_{i,A}}=37 \sim dgamma(\theta_B, 39,28)$

Ada dugaan bahwa banyaknya kasus korupsi di lembaga B lebih banyak dibandingkan dengan di lembaga A. Secara matematis, dugaan ini dapat dituliskan sebagai

$Pr(\theta_A<\theta_B|\sum{Y_{i,A}=36}, \sum{Y_{i,B}=37})=p$ ,

untuk suatu nilai $0\leq p\leq 1$ . Bagaimana menghitung nilai ? Dengan menyelesaikan integral berikut, $Pr(\theta_A<\theta_B|y_{1,A}, ..., y_{n_B,B})$ $=\int_0^{\infty}\int_0^{\theta_B}p(\theta_A, \theta_B|y_{1,A},...,y_{n_B,B} )d\theta_Ad\theta_B$

$=\int_0^{\infty}\int_0^{\theta_B}dgamma(\theta_A, 38,43)\times dgamma(\theta_B, 39,28)d\theta_A d\theta_B$

$=\frac{43^{38}28^{39}}{\Gamma(38)\Gamma(39)}\int_0^{\infty}\int_0^{\theta_B}\theta_A^{37}\theta_B^{38}e^{-43\theta_A-28\theta_B}d\theta_A d\theta_B$

Bagi yang tidak terlalu menyukai kalkulus , aproksimasi Monte Carlo bisa menjadi alternatif pilihan untuk mencari nilai dari integral di atas. Bagaimana caranya?

Misalkan kita tertarik terhadap parameter $\theta$ .
Misalkan pula adalah realisasi dari sampel yang berasal dari distribusi $p(y_1, ..., y_n|\theta)$ .
Seandainya kita bisa mengambil sampel, misal sebanyak nilai $\theta$ dari distribusi posterior $p(\theta|y_a, ..., y_n)$ , yaitu:

$\theta^{(1)}, ..., \theta^{(M)} \sim i.i.d. p(\theta|y_1, ..., y_n)$ Maka, distribusi empiris dari sampel $\left\{\theta^{(1)}, ..., \theta^{(M)} \right\}$ merupakan aproksimasi untuk $p(\theta|y_1, ..., y_n)$ , dengan tingkat 'keakuratan' makin tinggi untuk

yang lebih besar.

Law of large numbers

Jika $\theta^{(1)}, ..., \theta^{(M)$ adalah sampel i.i.d dari distribusi $p(\theta|y_1, ..., y_n)$ , maka $\frac{1}{M}\sum_{m=1}^Mg(\theta^{(m)})\rightarrow E[g(\theta|y_1, ..., y_n)]=\int g(\theta)p(\theta|y_1, ..., y_n)d\theta$ untuk $M\rightarrow \infty$ .

Yang berarti:

$\bar{\theta}=\sum_{m=1}^M\frac{\theta^{(m)}}{M}\rightarrow E[\theta|y_1,...,y_n]$
$\sum_{m=1}^M\frac{(\theta^{(m)}-\bar{\theta})^2}{M-1}\rightarrow Var[\theta|y_1, ..., y_n]$
$\frac{\text{banyaknya }(\theta^{(m)}\leq {c})}{M}\rightarrow Pr(\theta \leq c|y_1, ..., y_n)$
distribusi empiris dari $\left\{\theta^{(1)}, ..., \theta^{(M)}\right\} \rightarrow p(\theta|y_1, ..., y_n)$
median dari $\left\{\theta^{(1)}, ..., \theta^{(M)}\right\} \rightarrow \theta_{1/2}$

Dan berbagai bentuk informasi lain tentang $\theta$ dalam bentuk $g(\theta)$ yang berbeda, dapat diaproksimasi menggunakan sampel yang diperoleh menggunakan Monte Carlo. Misal, untuk distribusi posterior dari $\theta_A$ pada contoh di atas,

$\theta_A\sim dgamma(\theta_A, 38,43)$ .

Jika diambil sampel untuk $\theta_A$ menggunakan metode Monte Carlo, maka seperti dalam grafik berikut: [caption id="attachment_2439" align="aligncenter" width="300"]

Histogram dan densitas (empiris = merah, sebenarnya = biru) untuk data dari distribusi Gamma (38,43) menggunakan metode Monte Carlo. M adalah ukuran sampel.[/caption] Terlihat bahwa semakin besar ukuran sampel, maka distribusi empiris makin mendekati distribusi yang sebenarnya. Untuk contoh model Poisson di atas, karena dapat diperoleh solusi analitik untuk distribusi posteriornya:

$\theta_A\sim dgamma(\theta_A, 38,43)$ , dan

$\theta_B \sim dgamma(\theta_B, 39,28)$

maka, beberapa ringkasan numerik untuk kedua parameter tersebut juga dapat diperoleh langsung. Misal untuk data A:

mean $\theta_A$ $=\frac{38}{43}=0,8837$
$Pr(\theta_A<0.75)=\int_0^{0.75}dgamma(38,43)d\theta$
Interval kepercayaan $95\%$ sedemikian sehingga $\int_{LB}^{UB}dgamma(38,43)d\theta=0.95$

Menggunakan R,

Perbandingan ketiga nilai kuantitas di atas secara eksak dengan metode Monte Carlo dapat dirangkum dalam dalam tabel berikut:

Terlihat bahwa semakin besar ukuran sampel, maka taksiran Monte Carlo makin mendekati nilai yang sesungguhnya. Sebesar apakah ukuran sample yang diperlukan? Tergantung pada seberapa besar tingkat toleransi kita terhadap penyimpangan nilai taksiran terhadap nilai sesungguhnya. Bagaimana mengukurnya? Gunakan standar error Monte Carlo.

Central limit theorem:

Misalkan $\bar{\theta}=\sum_{m=1}^M\frac{\theta^{(m)}}{M}$ adalah mean dari sampel Monte Carlo, maka $\bar{\theta}\sim N (mean(\bar{\theta})=E[\theta|y_1, ..., y_n], Var(\bar{\theta})=\frac{Var[\theta|y_1, ... ,y_n]}{M})$

Sehingga, standar error Monte Carlo adalah

$\sqrt{\frac{\hat{\sigma}^2}{M}}$ , dengan $\hat{\sigma}^2=\sum\frac{(\theta^{(m)}-\bar{\theta})^2}{M-1}$ .

Aproksimasi Monte Carlo untuk interval kepercayaan 95% dari mean posterior $\theta$ adalah $\hat{\theta}\pm 2\sqrt{\frac{\hat{\sigma}^2}{M}}$ . Bagaimana menggunakannya? Misal, dari sampel Monte Carlo berukuran M=100 diperoleh $Var[\theta|y_1, ..., y_n]=0.36$ . Maka, standard error Monte Carlo adalah $\sqrt{0.36/100}=0.06$ . Jika ingin selisih antara hasil estimasi dari Monte Carlo dengan $E[\theta|y_1, ..., y_n]$ sekecil mungkin, misal <0.05

dengan probabilitas tinggi, maka ukuran sampel Monte Carlo harus ditambah menjadi:

$2\sqrt{\frac{.36}{M}}<0.05 \rightarrow M>676$ .

Secara grafis juga bisa ditentukan kira - kira ukuran sampel yang diperlukan supaya estimasi Monte Carlo sedekat mungkin dengan nilai sebenarnya. Dengan syntax R berikut,

diperolah grafik perubahan nilai mean, $Pr(\theta<1.75)$ , dan kuantil ke 97.5% ( $\theta_{\alpha}$ dimana $Pr(\theta<\theta_{\alpha})=0.975$ ) untuk ukuran sampel Monte Carlo hingga 1000. [caption id="attachment_2443" align="aligncenter" width="300"]

Nilai estimasi berdasarkan sampel Monte Carlo. Saat ukuran sampel membesar, estimasi semakin stabil. Garis horizontal adalah nilai eksak.[/caption] Bagaimana jika ingin melakukan prediksi? Gunakan distribusi predictive.

Kembali pada contoh model Poisson:

Jika diketahui nilai $\theta_A$ sebenarnya, maka dapat ditentukan probabilitas terjadi $\tilde{Y}$ kasus korupsi di lembaga A menggunakan :

Model sampling : $Pr(\tilde{Y}=\tilde{y}|\theta_A)=\frac{e^{-\theta_A}\theta_A^{\tilde{y}}}{\tilde{y}!}$ .

Akan tetapi, nilai kita tidak tahu nilai sebenarnya dari $\theta_A$ , sehingga diduga melalui prior $p(\theta)$ .

Untuk prediksi, selanjutnya berdasarkan apakah ada tambahan informasi data sampel atau tidak.

Seandainya tidak ada data sample, maka distribusi predictivenya adalah:

$Pr(\tilde{Y}=\tilde{y})=\int p(\tilde{y}|\theta_A)p(\theta_A)d\theta_A$

Bentuk di atas disebut prior predictive distribution.

Seandainya ada data sample $Y_1, ..., Y_{n_A}$ , maka distribusi predictive adalah:

$Pr(\tilde{Y}=\tilde{y}|Y_1=y_1, ..., Y_{n_A}=y_{n_A})=\int p(\tilde{y}|\theta_A, y_1, ...., y_{n_A})p(\theta_A|y_1, ...., y_{n_A})d\theta_A$

$=\int p(\tilde{y}|\theta_A)p(\theta_A|y_1, ...., y_{n_A})d\theta_A$

Bentuk di atas disebut posterior predictive distribution.

Dalam aplikasinya, $p(\tilde{y}|y_1, ..., y_n)$ cukup rumit untuk dievaluasi secara langsung. Prosedur Monte Carlo bisa digunakan dalam hal ini, secara tidak langsung.

$\text { ambil sampel }\theta^{(1)}\sim p(\theta|y_1, ..., y_n), \rightarrow \text{ ambil sampel } \tilde{y}^{(1)} \sim p(\tilde{y}|\theta^{(1)})$

$\text { ambil sampel }\theta^{(2)}\sim p(\theta|y_1, ..., y_n), \rightarrow \text{ ambil sampel } \tilde{y}^{(2)} \sim p(\tilde{y}|\theta^{(2)})$

...

$\text { ambil sampel }\theta^{(M)}\sim p(\theta|y_1, ..., y_n), \rightarrow \text{ ambil sampel } \tilde{y}^{(M)} \sim p(\tilde{y}|\theta^{(M)})$

$\left\{\tilde{y}^{(1)}, ..., \tilde{y}^{(M)}\right\}$ adalah sampel dari distribusi marginal posterior $\tilde{Y}$ , a.k.a distribusi posterior predictive. Penghitungan statistik yang akan diprediksi dengan mudah dilakukan pada data sampel ini. Sebagai contoh, kembali pada pertanyaan di awal pembahasan: Ada dugaan bahwa banyaknya kasus korupsi di lembaga B lebih banyak dibandingkan dengan di lembaga A. Secara matematis, dugaan ini dapat dituliskan sebagai $Pr(\theta_A<\theta_B|\sum{Y_{i,A}=36}, \sum{Y_{i,B}=37})=p$ , dengan suatu nilai $0\leq p\leq 1$ . Bagaimana menghitung nilai ?

$\theta^{(1)}\sim gamma(a+\sum y_i , b+n), \quad \rightarrow \quad \tilde{y}^{(1)} \sim Poisson(\theta^{(1)})$

$\theta^{(2)}\sim gamma(a+\sum y_i , b+n), \quad \rightarrow \quad \tilde{y}^{(2)} \sim Poisson(\theta^{(2)})$

...

$\theta^{(M)}\sim gamma(a+\sum y_i , b+n), \quad \rightarrow \quad \tilde{y}^{(M)} \sim Poisson(\theta^{(M)})$

Dengan menggunakan R, prosedur di atas adalah:

Hasil di atas menunjukkan bahwa probabilitas mendapatkan lebih banyak kasus korupsi di lembaga B dibandingkan dengan lembaga A adalah 0.3457?

Apakah angka ini cukup menunjang klaim kebijakan e-budgeting menekan angka korupsi?

Syntax R untuk contoh di atas dapat diunduh di sini(montecarlo).

Aproksimasi Monte Carlo

No comments:

Post a Comment

SEARCH

LATEST

FOLLOW ME

Visitors

Translate~Terjemahkan

Pages

Follow Us

Popular

Archive

Postingan Populer

Latest courses

Comments

About

Top Links Menu