Inferensi Bayesian menggunakan Model Poisson ~ tau-data Indonesia

Misal peubah acak berdistribusi Poisson $(\theta)$ . Fungsi densitasnya diberikan oleh $Pr(Y=y|\theta)=dpois(y,\theta)=\theta^y \frac{e^{-\theta}}{y!}$ untuk $y\in \{0,1,2,...\}$

dengan

mean $E[Y|\theta]=\theta$ dan variansi $Var [Y|\theta]}=\theta$ .

Sehingga pada data yang berdistribusi Poisson, keragaman data akan membesar seiring dengan kenaikan mean data.

Bagaimana melakukan inferensi secara Bayesian menggunakan model Poisson?

Prinsip dasar metode Bayesian secara umum adalah:

Sehingga, untuk melakukan inferensi, diperlukan asumsi awal (prior) dan tambahan informasi dari data (likelihood) untuk kemudian menghasilkan informasi terbaharukan (update).

Maka, selanjutnya kita bahas dahulu komponen – komponen tersebut.

Likelihood:

Untuk sampel berukuran , jika $Y_1,Y_2, ..., Y_n \sim i.i.d\ Poisson(\theta)$ , maka fungsi likelihoodnya adalah:

$Pr(Y_1=y_1, ..., Y_n=y_n|\theta)=\prod_{i=1}^np(y_i |\theta) =\prod_{i=1}^n\frac{1}{y_i!}\theta ^{y_i}e^{-\theta} =k(y_1, ...., y_n)\theta^{\sum{y_i}}e^{-n\theta}$ (1)

Prior
Salah satu kelas prior yang memudahkan dalam penghitungan posterior adalah prior conjugate , karena posterior maupun prior berasal dari kelas distribusi yang sama. Untuk pembahasan kali ini akan digunakan conjugate prior.

Bagaimana menentukan prior conjugate ?

Posterior dari $\theta$ dengan prior $p(\theta)$ adalah:

$Posterior\propto {\color{Red} prior} \times {\color{Blue} likelihood}\\ p(\theta|y_1, ..., y_n) \propto{\color{Red} p(\theta)}\times {\color{Blue} p(y_1, ..., y_n|\theta)}\propto {\color{Red} p(\theta)}\times {\color{Blue} \theta^{\sum{y_i}}e^{-n\theta}}$

Dari bentuk di atas, karena distribusi prior conjugate berasal dari kelas yang sama dengan distribusi posterior, maka bentuk distribusi prior haruslah memuat komponen $\theta^{k_1}e^{-k_2\theta}$ untuk suatu konstanta dan .

Salah satu kelas distribusi yang memenuhi bentuk tersebut adalah distribusi Gamma.

Sehingga, pilih distribusi prior untuk $\theta$ :

$\theta\sim dgamma (a,b)$ ,

dengan fungsi densitas

$p(\theta)=\frac{b^a}{\Gamma(a)}\theta^{{\color{Red} a}-1}e^{{\color{Red} -b}\theta} \,\text{untuk}\ \theta, a, b>0$ (2)

Mean dan variansi dari distribusi Gamma diberikan oleh

$E[\theta]=\frac{a}{b}$ dan $Var[\theta]=\frac{a}{b^2}$ .

Beberapa bentuk densitas distribusi Gamma ditampilkan berikut:

Posterior

Dari pembahasan sebelumnya, dengan

Data: $Y_1, ..., Y_n |\theta \sim Poisson(\theta)$ , i.i.d, dengan likelihood seperti (1), dan

Prior untuk $\theta$ : $p(\theta)$ dengan densitas seperti pada (2), maka

Posterior untuk $\theta$ :

$p(\theta|y_1, ..., y_n)= p(\theta)\times \frac{p(y_1, ..., y_n|\theta)}{p(y_1, ..., y_n )} \\ =\theta^{a-1} e^{-b\theta}\times\theta^{\sum{y_i}} e^(-n\theta)\times k(y_1, ...,y_n,a,b)\\ =[\theta^{{\color{Red} a+\sum{y_i}}-1} e^{{\color{Red} -(b+n)}\theta)} ]\times [k(y_1, ..., y_n,a,b)]=dgamma({a+\sum_{i=1}^n {Y_i}}, b+n)$ Sehingga, untuk model Poisson dengan:
Prior: $\theta \sim dgamma \,(a,b)$ ,

Likelihood data: $Y_1, ..., Y_n |\theta \sim dPoisson (\theta)$ ,

maka

Posterior:

$\theta|Y_1, ..., Y_n \sim dgamma (a+\sum_{i=1}^n Y_i, b+n)$

Mean posterior dari $\theta$ adalah:

$E[\theta|y_1, ..., y_n ]=\frac{a+\sum{y_i}}{b+n}=\frac{b}{b+n}\times\frac{a}{b}+\frac{n}{b+n}\frac{\sum{y_i}}{n}$

merupakan jumlah terboboti antara mean prior $\frac{a}{b}$ dan mean sampel $\frac{\sum{y_i}}{n}$ .

Analogi dengan ukuran sampel dan jumlah nilai data $\sum{y_i}$ , maka

dapat diinterpretasikan sebagai banyaknya pengamatan prior, dan
sebagai penjumlahan nilai dari pengamatan prior tersebut.

Jika ukuran sampel , maka

$E[\theta|y_1, ..., y_n ]\approx \frac{\sum{y_i}}{n}=\bar{y}$

$Var[\theta|y_1, ..., y_n]\approx \frac{\bar{y}}{n}$ yang berarti informasi posterior akan didominasi oleh infomasi dari data sampel.

Bagaimana melakukan prediksi untuk data baru?

Gunakan distribusi posterior predictive.

Secara matematis, distribusi posterior predictive dinyatakan oleh $p(\tilde{y}|y_1, ..., y_n )$ dengan $\tilde{y}$ adalah data baru.

Distribusi posterior predictive untuk model di atas adalah:

$p(\tilde{y}|y_1,...,y_n )=\int_\theta p(\tilde{y}|\theta,y_1, ....,y_n )p(\theta|y_1, ....,y_n )d\theta〗=∫p(\tilde{y}|\theta)p(\theta|y_1, ...., y_n )d\theta$

$=\int dpois(\tilde{y},\theta)dgamma(\theta,a+\sum{y_i},b+n)d\theta$

$=\int\left \{ \frac{1}{\tilde{y!}}\theta^{\tilde{y}}e^{-\theta} \right\}\left\{\frac{(b+n)^{a+\sum{y_i}}}{\Gamma(a+\sum{y_i)}} \theta^{a+\sum{y_i}-1}e^{-(b+n)\theta}\right \}d\theta$

$=\frac{(b+n)^{a+\sum{y_i}}}{\Gamma(\tilde{y}+1)\Gamma(a+\sum{y_i})}\int_{\theta}\theta^{a+\sum{y_i}+\tilde{y}-1}e^{-(b+n+1)\theta}d\theta$

Dengan menggunakan

$1=\int_0^{\infty}\theta^{a-1}e^{-b\theta}d\theta=\frac{\Gamma(a)}{b^a}$ untuk

Maka:

$p(\tilde{y}|y_1,...,y_n )=\frac{\Gamma(a+\sum{y_i}+\tilde{y} )}{\Gamma(\tilde{y}+1)\Gamma(a+\sum{y_i})}(\frac{b+n}{b+n+1})^{a+\sum{y_i}} (\frac{a}{b+n+1})^{\tilde{y}}$

$\sim NB(a+\sum{y_i}, b+n)$ untuk $\tilde{y} \in \left\{ 0, 1, 2, ...\right\}$ .

Contoh

Data mengenai banyaknya kasus korupsi di dua institusi pemerintah, sebut lembaga A dan lembaga B. Lembaga A menerapkan prinsip e-budgeting, sementara lembaga B tidak. Dari 40 divisi pada lembaga A ditemukan 36 kasus korupsi, sedangkan dari 25 divisi pada lembaga B, ditemukan 37 kasus serupa. Ringkasan data pada grafik berikut.

Misal akan dibandingkan kedua kementrian tersebut berdasarkan jumlah kasus korupsi yang berhasil diungkap.

Maka, modelnya adalah:

$Y_{1,A}, ..., Y_{n_A,A} |\theta_A \sim i.i.d. Poisson (\theta_A )$ , untuk lembaga A

$Y_{1,B}, ..., Y_{n_B,B} |\theta_B \sim i.i.d. Poisson (\theta_B )$ , untuk lembaga B.

Ringkasan numerik untuk kedua data tersebut:

$n_A=40,\quad \sum{Y_{i,A}}=36,\quad \bar{Y_A}=0.9$

$n_B=25,\quad \sum{Y_{i,B}}=37,\quad \bar{Y_B}=1.48$

Jika parameter $\left\{\theta_A, \theta_B\right\}\sim i.i.d \,dgamma(a=2, b=3)$ , maka distribusi posterior dari kedua parameter tersebut adalah:

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

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

Taksiran posterior untuk mean, modus, dan interval kepercayaan 95% untuk dan dapat diperoleh dari distribusi Gamma posterior.

Grafik distribusi prior (sama untuk A dan B) dan posterior dapat diperoleh sebagai berikut:

Dan diperoleh hasil

Distribusi dari parameter theta (hitam = prior, merah dan biru = posterior untuk A dan B)

Intepretasi:

Terlihat bahwa likelihood cukup berperan meng-update informasi dari prior, sebagaimana ditunjukkan oleh distribusi posterior (A = merah, B = biru) yang berbeda dengan distribusi prior (hitam). Grafik tersebut juga mengindikasikan bahwa jumlah kasus korupsi secara rata – rata ( di lembaga A relative lebih sedikit daripada di lembaga B (grafik merah di sebelah kiri grafik biru ().

Tapi itu kan secara rata – rata, gambaran umum. Bagaimana untuk prediksi, kasus per kasus? Apakah lembaga A , dengan kebijakan e-budgeting, memang memiliki kasus korupsi yang lebih sedikit dibanding lembaga B?

Untuk menjawab hal ini, kita dapat menghitung peluang terjadinya hal tersebut menggunakan posterior predictive distribution. Sebagaimana telah ditunjukkan di atas, distribusi binomial negative merupakan posterior predictive distribution dari $\tilde{Y_A}$ dan $\tilde{Y_B}$ .

Misalkan akan dilihat berapa peluang terjadi kasus korupsi, untuk .

pA adalah peluang terjadinya kasus korupsi di lembaga A, dan pB untuk hal yang sama di lembaga B. Penyajian grafis di bawah ini:

picPois10 — Posterior predictive probabilitas untuk ditemukannya kasus korupsi pada lembaga A (hitam) dan lembaga B (merah).

Intepretasi:

Dari plot di atas terlihat bahwaprobabilitas untuk tidak ditemukannya kasus korupsi di lembaga A jauh lebih tinggi daripada lembaga B, hampir sama untuk 1 kasus korupsi, dan untuk lebih dari 1 kasus lebih berpeluang untuk terjadi di lembaga B.

Data untuk contoh dapat diunduh di sini (data-korupsi).

Syntax R dapat diunduh di sini (model-poisson_syntax).

Share this:

Related