Secara umum, masyarakat ilmiah belum sepenuhnya bisa memanfaatkan teknik Bayesian untuk melakukan inferensi ketika dihadapkan pada masalah nyata. Hal ini terutama disebabkan oleh dua permasalahan yang berbeda, yaitu ketidakmampuan untuk menghitung integral dimensi tinggi yang diperlukan untuk mengkarakterisasi distribusi posterior pada model. Sebagian besar, hal ini dapat diselesaikan dengan munculnya Markov Chain Monte Carlo (MCMC).
Metode MCMC secara luas dianggap sebagai perkembangan yang paling penting dalam komputasi statistik dalam sejarah. MCMC memungkinkan ahli statistik untuk mengepas beberapa model dasar probabilitas termasuk yang bahkan tidak dianggap beberapa tahun yang lalu. Sayangnya, statistik satu-satunya yang telah bersedia dan mampu menulis kode komputer yang diperlukan untuk menggunakan metode MCMC sesuai model probabilitas.
Salah satu kegunaan dari teori rantai Markov ergodik adalah Markov Chain Monte Carlo atau MCMC. Misalkan akan dilakukan simulasi dari suatu probabilitas (atau disebut juga densitas target)permasalahan muncul karena simulasi secara langsung tidak mungkin dilakukan atau infeasible dikarenakan dimensi yang tinggi. Maka MCMC membantu menyelesaikan masalah ini, Permasalahan semacam ini seringkali muncul dalam berbagai macam aplikasi sains contohnya pada bidang statistika, komputer sains, dan fisika statistik. Markov Chain Monte Carlo memberikan solusi secara tidak langsung pada permasalahan yang didasarkan pada observasi dengan membangun sebuah rantai Markov ergodik dan sebagai ukuran probabilitas stasioner. Ini lebih efektif dibandingkan melakukan simulasi secara langsung pada (Robert & Rosenthal,2003).
Sebuah barisan $X_1,X_2,\cdots ,X_n$ elemen acak dari beberapa himpunan adalah Markov Chain jika distribusi kondisional $X_{n+1}$ diberikan $X_1, ,\cdots ,X_n$ bergantung hanya pada $X_n$. Dimana Himpunan bernilai $X_i$ disebut state space dari Markov Chain.
Sebuah Markov Chain mempunyai stasionary transition probabilities jika distribusi kondisional dari $X_{n+1}$ diberikan $X_{n}$ tidak bergantung pada $n$. Ini adalah jenis utama dari Markov Chain yang menarik di dalam MCMC. Distribusi bersama dari markov chain ditentukan oleh :
1. Distribusi marginal dari $X_i$ disebut Initial Distribution.
2. Distribusi kondisional dari $X_{n+1}$ diberikan $X_{n}$, disebut transition probability distribution (karena asumsi stationary transition probabilities ini tidak bergantung pada $n$).
Seseorang memperkenalkan Markov Chain melalui tipe khusus pada proses stokastik yang biasanya hanya melihat contoh dimana state space nya terbatas atau bisa dihitung. Jika state space terbatas, maka dapat dituliskan $\left \{ X_1,\cdots ,X_n \right \}$, maka initial distribution dapat dikaitkan dengan sebuah vector $\lambda =\left \{ \lambda _1,\cdots ,\lambda _n \right \}$ yang dapat didefinisikan dengan
Dan transition distribution dapat dikaitkan dengan sebuah matriks $P $mempunyai elemen $P_{ij}$ didefinisikan oleh
Ketika keadaan ruang adalah dapat dihitung tak terbatas, kita bisa memikirkan sebuah vektor yang tak terbatas dan matriks. Tapi kebanyakan rantai Markov yang menarik di MCMC memiliki keadaan ruang tak terhitung, dan kemudian kita tidak bisa memikirkan distribusi awal sebagai vektor atau distribusi probabilitas transisi sebagai sebuah matriks. Kita harus menganggap mereka sebagai distribusi probabilitas bersyarat dan distribusi probabilitas bersyarat.
Algoritma Monte Carlo adalah metode Monte Carlo numerik yang digunakan untuk menemukan solusi problem matematis (yang dapat terdiri dari banyak variabel)yang susah dipecahkan, misalnya dengan kalkulus integral, atau metode numerik lainnya. Pada tahun 1951 Von Neumann memberikan suatu metode untuk membangkitkan sampel dari sebuah densitas $f=\left (\frac{1}{K} \right )f_{1}$ hanya diketahui $f_{1}$ jika ada suatu densitas $h$ sedemikian sehingga (adalah mungkin untuk membangkitkan sampel dari $h$ dan) $ f_{1}(x)\leq Mh(x)$ $ \forall_{x}$ Algoritma diberikan di bawah ini yang dikenal sebagai rejection sampling :
Bangkitkan sampel acak dari distribusi dengan densitas $h$ Misal $y$.
Terima $y$ sebagai sampel dengan probabilitas $ \left [\frac{f_1(y)}{Mh(y)} \right ]$.
Jika langkah ke 2 tidak berhasil kembali ke langkah 1.
Markov Chain Monte Carlo (MCMC)adalah sebuah metode untuk membangkitkan peubah-peubah acak yang didasarkan pada rantai markov. Dengan MCMC akan diperoleh sebuah barisan sampel acak yang berkorelasi, yakni nilai ke–j dari barisan $ \left \{ \theta _j \right \}$ disampling dari sebuah distribusi peluang yang bergantung pada nilai sebelumnya $ \left \{ \theta _{j-1} \right \}$. Distribusi eksak dari $ \left \{ \theta _j \right \}$ umumnya tidak diketahui, namun distribusi pada setiap iterasi dalam barisan nilai sampel tersebut akan konvergen pada distribusi yang sesungguhnya untuk nilai j yang cukup besar. Oleh karena itu, jika ukuran sampel yang diperbaharui cukup besar maka kelompok terakhir dari nilai yang disampling dalam barisan tersebut, misal $\left \{ \theta _{p+1},\theta _{p+2},\cdots \right \}$ akan mendekati sebuah sampel yang berasal dari distribusi yang diinginkan. Notasi P biasanya disebut sebagai burn in period.
Terdapat dua algoritma utama dalam MCMC, yaitu algoritma Metropolis-Hastings dan algoritma Gibbs Sampling.
Algoritma Metropolis-Hasting (MH)
Algoritma Metropolis-Hasting (MH)digunakan untuk membantu membangkitkan sampel-sampel acak dari distribusi posterior yang diinginkan. Dalam algoritma MH dibutuhkan sebuah proposal distribution $p(\theta |\theta _{j-1})$ untuk membangkitkan kandidat sampel acak. Langkah-langkah dasar dari algoritma ini adalah
Langkah 1 : Ambil nilai awal, yaitu $\theta _0$. Untuk iterasi j = 1. Bangkitkan $\theta ^{*}\sim p(\theta |\theta _{j-1})$.
Langkah 2 : Bangkitkan sampel acak u dari distribusi uniform $ U [ 0 ,1 ]$.
Langkah 3 : Jika $ u<min\left (1, \frac{p\left (\theta ^*|X,y \right )p\left (\theta _{j-1}|\theta ^* \right )}{p\left ( \theta_{j-1}|X,y \right )p\left (\theta^*|\theta_{j-1} \right )}\right ),$ ambil $\theta_j=\theta ^*$ . Namun jika $ u<min\left (1, \frac{p\left (\theta ^*|X,y \right )p\left (\theta _{j-1}|\theta ^* \right )}{p\left (\theta_{j-1}|X,y \right )p\left (\theta^*|\theta_{j-1} \right )}\right ),$ ambil $\theta_j=\theta _{j-1}$
Langkah 4 : Ulangi langkah 1-3 sampai dengan jumlah yang diinginkan
Statistik yang digunakan untuk mengukur derajat ketergantungan antar pengambilan secara berurutan dalam sebuah rantai markov adalah autokorelasi. Autokorelasi mengukur korelasi antara dua buah himpunan dari nilai-nilai yang tersimulasi $\left \{ \theta _j \right \}$ dan $\left \{ \theta _{j+L} \right \}$ dengan L adalah lag atau bilangan yang memisahkan dua himpunan tersebut. Untuk sebuah hyperparameter $\theta _j$ tertentu, nilai autokorelasi pada lag ke-L dapat dihitung dengan formula
dengan M adalah ukuran sampel acaknya.
Algoritma Gibbs Sampling
Gibbs Sampling bisa diterapkan apabila distribusi probabilitas bersama (joint probability distribution)tidak diketahui secara eksplisit, tetapi distribusi bersyarat (conditional distribution)dari tiap-tiap variabel diketahui. Algoritma Gibbs sampling bisa dituliskan sebagai berikut :
Langkah 1 : Tentukan nilai awal $ x^0=\left (x _{1}^{\left (0 \right )},\cdots , x _{p}^{\left (0 \right )}\right )$
Langkah 2 : Ulangi Langkah 2 untuk $j=1,2,\cdots ,M$
Bangkitkan $X_{1}^{(j+1)}$ dari $f_1\left (x_1|x_{2}^{\left (j \right )},x_{3}^{\left (j \right )},\cdots ,x_{p}^{\left (j \right )} \right )$
Bangkitkan $X_{2}^{(j+1)}$ dari $f_2\left (x_2|x_{1}^{\left (j +1\right )},x_{3}^{\left (j \right )},\cdots ,x_{p}^{\left (j \right )} \right )$
Bangkitkan $X_{p}^{(j+1)}$ dari $f_p\left ( x_p|x_{1}^{\left (j +1\right )},x_{3}^{\left (j \right )},\cdots ,x_{p-1}^{\left (j \right )} \right )$
Langkah 3 : Kembalikan nilai $\left \{ x^1, x^2,\cdots ,x^M \right \}$
densitas $f_1, f_2,\cdots ,f_p$ disebut distribusi bersyarat penuh, dan densitas yang digunakan untuk simulasi. Walaupun dalam dimensi tinggi semua simulasi adalah univariate. Dalam gibbs sampling tidak ada mekanisme penerimaan dan penolakan semua sampel hasil simulasi diterima.
MCMCregress membangkitkan sampel dari distribusi posterior dari model regresi linier dengan error berdistribusi Gaussian menggunakan Gibbs sampling (dengan prior multivariat gaussian pada vektor beta, dan prior Invers Gamma pada varians error). Persediaan pengguna data dan prior, dan sampel dari distribusi posterior dikembalikan sebagai objek MCMC, yang dapat kemudian dianalisis dengan fungsi yang disediakan dalam paket coda.
Dengan bentuk model yang digunakan adalah
$y_i=x_i'\beta + \epsilon_i,$
dengan error diasumsikan berdistribusi Gaussian $N(0,\sigma^2).$
Asumsi standar prior semi conjugate $\beta\sim N(b0,B0^{-1})$ dan $\sigma^{-2} \sim \gamma(\frac{c0}{2},\frac{d0}{2})$.
dengan $\beta$ dan $\sigma^{2}$ diasumsikan prior independent.
Syntax R untuk Regresi MCMC pada data Gaussian adalah sebagai berikut :
MCMCregress(formula, data = NULL, burnin = 1000,
mcmc = 10000, thin = 1, verbose = 0, seed = NA,
beta.start = NA, b0 = 0, B0 = 0, c0 = 0.001,
d0 = 0.001, sigma.mu = NA, sigma.var = NA,
marginal.likelihood = c("none", "Laplace", "Chib95"), ...)
Fungsi ini membangkitkan sampel dari distribusi posterior dari model regresi Poisson menggunakan algoritma random walk Metropolis. informasi prior dan sampel dari distribusi posterior dikembalikan sebagai objek MCMC, yang kemudian dapat dianalisis dengan fungsi yang disediakan dalam paket. Dengan bentuk model yang digunakan adalah $y_i \sim Poisson(\mu _i)$ Dimana fungsi link yang digunakan adalah $\mu _i = \exp(x_i' \beta )$. Sedangkan distribusi prior-nya diasumsikan berdistribusi beta, yaitu $\beta\sim Be(b0,B0^{-1})$ .
Syntax R untuk Regresi MCMC pada data Poisson adalah sebagai berikut :
MCMCpoisson(formula, data = NULL, burnin = 1000, mcmc = 10000,
thin = 1, tune = 1.1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, marginal.likelihood = c("none", "Laplace"),...)
MCMClogit mensimulasikan sampel dari distribusi posterior pada model regresi logistik menggunakan algoritma random walk Metropolis. Model yang digunakan adalah $ y_i\sim Bernoulli(\pi _i)$ Dimana fungsi link yang digunakan adalah $ \pi _i =\frac{\exp(x_i' \beta )} {[1 + \exp(x_i'\beta )]}$ Sedangkan distribusi prior-nya diasumsikan berdistribusi beta, yaitu $\beta \sim Be(b0,B0^{-1})$. Syntax R untuk Regresi MCMC pada distribusi logistic adalah sebagai berikut :
MCMClogit(formula, data=NULL, burnin = 1000, mcmc = 10000,
thin=1, tune=1.1, verbose = 0, seed = NA, beta.start = NA,
b0 = 0, B0 = 0, user.prior.density=NULL, logfun=TRUE,
marginal.likelihood=c("none", "Laplace"), ...)
Analisis Bayesian menggali dua buah sumber informasi tentang parameter suatu model statistik. Sumber informasi pertama berasal dari sampel yang disebut informasi sampel (sampel information ). Sumber informasi kedua berasal dari opini sang ahli yang disebut informasi prior (prior information ). Gabungan dua buah sumber informasi ini akan membentuk informasi posterior (posterior information ). Penggabungan kedua sumber informasi ini dicapai melalui teorema Bayes (Bayes’ theorem ). Misalkan parameter yang dimaksud adalah vektor $\theta$ dan sampelnya adalah vektor $y$. Parameter $\theta$ dan sampel $y$. keduanya adalah variabel acak (random variables)dengan fungsi densitas gabungan $h\left (\theta ,y \right )$. Observasi sampel y berasal dari fungsi densitas $f\left (y| \theta \right )$ yang identik dengan fungsi likelihood $l\left (\theta |y \right )$. Fungsi likelihood $l\left (\theta |y \right )$ ini memuat seluruh informasi sampel tentang $\theta$. Fungsi densitas gabungan $h\left (\theta ,y \right )$ dapat dinyatakan dengan $$h\left (\theta|y\right )=g \left (\theta|y\right )f\left (y \right )=f\left (y| \theta\right )g\left (\theta \right )$$ Sehingga diperoleh $$g \left (\theta|y\right )= \frac{f\left (y| \theta\right )g\left (\theta \right )}{f\left (y \right )}$$ $$= \frac{l\left ( \theta|y\right )g\left (\theta \right )}{f\left (y \right )}$$ Persamaan (3)dikenal dengan teorema Bayes. Teorema Bayes ini lebih dikenal lewat ungkapan $g\left ( \theta|y\right )\propto l\left ( \theta|y\right )g\left (\theta \right )$ dengan notasi $\propto$ menyatakan ‘proporsional terhadap’. Hal ini dimungkinkan karena fungsi $f (y)$ tidak mengandung $\theta$ dan dapat dianggap sebagai konstanta. Fungsi $g (\theta |y)$ akan disebut fungsi densitas posterior (posterior density function)dan fungsi $g (\theta )$ akan disebut fungsi densitas prior (prior density function).
Distrribusi prior dikelompokkan menjadi dua berdasarkan bentuk fungsi likelihood, yaitu:
Berkaitan dengan bentuk distribusi hasil identifikasi pola datanya
Distribusi prior konjugat (conjugate), mengacu pada acuan analisis model terutama dalam pembentukan fungsi likelihoodnya sehingga dalam penentuan prior konjugat selalu dipikirkan mengenai penentuan pola distribusi prior yang mempunyai bentuk konjugat dengan fungsi densitas peluang pembangun likelihoodnya.
Distribusi prior tidak konjugat (non-conjugate), pemberian prior pada model tidak mempertimbangkan pola pembentuk fungsi likelihoodnya
Berkaitan dengan penentuan parameter pada pola distribusi prior
Distribusi prior informatif, mengacu pada pemberian parameter dari distribusi prior yang telah dipilih baik distribusi prior konjugat atau tidak, pemberian nilai parameter pada distribusi prior ini didasarkan pada informasi yang diperoleh
Distribusi prior non informatif, pemilihannya tidak didasarkan pada data yang ada atau distribusi prior yang tidak mengandung informasi tentang parameter θ. Apabila pengetahuan tentang priornya sangat lemah, maka bisa digunakan prior berdistribusi normal dengan mean nol dan varian besar. Efek dari penggunaan prior dengan mean nol adalah estimasi parameternya dihaluskan menuju nol. Pemulusan ini dilakukan oleh varian, sehingga pemulusan tersebut bisa dilakukan dengan meningkatkan varian
Distribusi posterior adalah fungsi densitas bersyarat θ jika diketahui nilai observasi x dan dapat ditulis sebagai berikut :
$f(\theta |y)=\frac{f(\theta ,y)}{f(y)}$
Fungsi kepadatan bersama dan marginal yang diperlukan dapat ditulis dalam bentuk distribusi prior dan fungsi likelihood,
$f(\theta ,y)=f(y|\theta )f(\theta )$
$f(y)=\int_{-\infty}^{\infty}f(\theta ,y)d\theta = \int_{-\infty}^{\infty}f(\theta )f(y|\theta )d\theta$
Sehingga fungsi densitas posterior untuk variabel random kontinu sebagai berikut,
$f(\theta |y)=\frac{f(\theta )f(y|\theta )}{\int_{-\infty}^{\infty}f(\theta )f(y|\theta )d\theta }$
Pada bagian ini disajikan secara naratif aktivasi data, langkah-langkah pemeriksaan asumsi, pemeriksaan GOF model sebagaimana teori yang disampaikan sebelumnya. Anda dapat mengaktifkan data dari database R, atau data yang anda miliki dalam format khusus (Tex atau CSV). Jika anda memiliki data dalam format excel, untuk saat ini anda harus mengkonversinya ke bentuk teks atau CSV terlebih dahulu.
| Pilihan Data |
Dari summary data yang ada, kita bisa menentukan variabel-variabel yang akan dijadikan variabel bebas (eksplanatori)dan variabel terikat (respon). Sebelum menentukan variabel bebas dan terikat (respon), kita dapat juga membuat matriks korelasi untuk mendeteksi variabel-variabel yang terindikasi memiliki hubungan. Hasil matriks korelasinya, atau matrks diagram pencarnya adalah sebagai berikut.
Luaran 1. Data Aktif Ukuran sampel untuk dianalisis dengan MCMC (<20)
Tampilan yang lebih rinci dari data dapat dilihat pada Lampiran.
Dari pilihan data yang dilakukan, kita bisa memperoleh gambaran kasar tentang hubungan $X$ dengan $Y$ (seperti kelinieran dan ada tidaknya pencilan)dengan melihat Grafik diagram pencarnya.
Untuk mendapatkan gambaran tentang distribusi $Y$ kita dapat memeriksa Grafik QQ-Plot, Histogram dan Box-Plot seperti berikut.
Jenis MCMC sesuai jenis data:
Keterangan :
"Gaussian" untuk jenis data kontinu, "Poisson" untuk jenis data berupa cacahan,
sedangkan "Binomial" untuk jenis data berupa cacahan 0 dan 1 (biner)
Banyaknya pemanasan/burned-in . Banyaknya MCMC .
Keterangan :
Pengujian hipotesis terhadap parameter regresi dilakukan dengan pendekatan interval konfidensi 95% dari masing-masing parameter. Hal ini dikarenakan distribusi posterior tidak diketahui dengan pasti. Interval konfidensi 95% dihitung dengan batas bawah yaitu kuantil ke 2,5% dan batas atasnya adalah kuantil ke 97,5%. Parameter dinyatakan signifikan jika interval konfidensi 95% parameter tidak memuat nilai nol.
MCMC adalah suatu metode simulasi yang me-monte-carlo-kan nilai parameter model yang sesuai dengan proses Markov Chain untuk mendapatkan data sampel berdasarkan skenario sampling tertentu.
Terdapat dua algoritma utama dalam MCMC, yaitu algoritma Metropolis-Hastings dan algoritma Gibbs Sampling.
Pada Bayesian semua parameter dalam model diperlakukan sebagai variabel.
Inferensi Bayesian dilakukan dengan menggunakan distribusi sampel dan distribusi awal (prior)yaitu dugaan distribusi populasi untuk mendapatkan distribusi poeterior dan estimasi parameter dilakukan pada distribusi posterior tersebut.