Berapa banyak sisi yang dimiliki die? Kesimpulan Bayesian dalam JAGS

Masalah

Saya ingin melakukan beberapa kesimpulan tentang sistem analog dengan mati dengan jumlah sisi yang tidak diketahui. Die digulung beberapa kali, setelah itu saya ingin menyimpulkan distribusi probabilitas atas parameter yang sesuai dengan jumlah sisi yang dimiliki die, θ.

Intuisi

Jika setelah 40 gulungan Anda telah mengamati 10 merah, 10 biru, 10 hijau dan 10 kuning, tampaknya θ harus memuncak pada 4, dan bias dari setiap sisi bergulir adalah distribusi yang berpusat pada 1/4.

θ memiliki batas bawah yang sepele, menjadi jumlah sisi yang berbeda yang diamati dalam data.

Batas atas masih belum diketahui. Mungkin ada sisi kelima yang mungkin memiliki bias yang rendah. Semakin banyak data yang Anda amati tanpa kategori kelima, semakin tinggi probabilitas posterior of = 4.

Pendekatan

Saya telah menggunakan JAGS untuk masalah yang sama (via R dan rjags) yang sepertinya cocok di sini.

Sehubungan dengan data, katakanlah obs <- c(10, 10, 10, 10)sesuai dengan pengamatan pada contoh di atas.

Saya pikir pengamatan harus dimodelkan dengan distribusi multinomial obs ~ dmulti(p, n), di mana p ~ ddirch(alpha)dan n <- length(obs).

θ ditautkan dengan jumlah kategori yang tersirat alpha, jadi bagaimana saya bisa memodelkan alphauntuk mencakup berbagai kategori yang berbeda?

Alternatif?

Saya cukup baru dalam analisis Bayesian, jadi mungkin menyalak sepenuhnya, apakah ada model alternatif yang mungkin memberikan wawasan berbeda tentang masalah ini?

Terimakasih banyak! David

Masalah yang menarik ini disebut 'pengambilan sampel spesies', yang telah menerima banyak perhatian selama bertahun-tahun, dan mencakup banyak masalah estimasi lainnya (seperti mark-recapture). Cukuplah untuk mengatakan, JAGS tidak akan membantu Anda dalam hal ini - JAGS tidak dapat menangani rantai Markov dengan dimensi variabel di seluruh iterasi. Seseorang harus beralih ke skema MCMC yang dirancang untuk masalah seperti itu, seperti MCMC lompat terbalik.

Inilah satu pendekatan yang cocok untuk model spesifik yang Anda gambarkan, yang pertama kali saya temui dalam karya Jeff Miller ( arxived ).

Bagian I (Pertanyaan awal)

Satu asumsi yang akan saya buat adalah bahwa pengamatan kategori tertentu menyiratkan adanya kategori dengan peringkat yang lebih rendah. Artinya, mengamati gulungan mati di sisi 9 menyiratkan keberadaan sisi 1-8. Tidak harus seperti ini - kategorinya bisa sewenang-wenang - tetapi saya akan berasumsi demikian dalam contoh saya. Ini berarti bahwa 0 nilai dapat diamati, berbeda dengan masalah estimasi spesies lainnya.

Katakanlah kita memiliki sampel multinomial,

$$Y = \{y_1, y_2, \dots, y_m, y_{m+1}, \dots, y_{n} \} \sim \mathcal{M}(\{p_1, p_2, \dots, p_m, p_{m+1}, \dots, p_n\})$$

$m$$n$$\{y_{m+1},\dots,y_{n}\}$$n$$[1, \infty)$

$$n \sim \mathcal{P}(\lambda), n > 0$$

Sebelum nyaman untuk probabilitas multinomial adalah Dirichlet,

$$P = \{ p_1, \dots, p_n \} \sim \mathcal{D}(\{ \alpha_1, \dots, \alpha_n \})$$

$\alpha_1 = \alpha_2 = \dots = \alpha_n = \tilde{\alpha}$

Untuk membuat masalah lebih mudah ditangani, kami meminggirkan bobot:

$$p(Y|\tilde{\alpha}, n) = \int_P p(Y|P, n)p(P|\tilde{\alpha}, n) dP$$

Yang dalam hal ini memimpin distribusi Dirichlet-multinomial yang dipelajari dengan baik . Tujuannya kemudian untuk memperkirakan posterior bersyarat,

$$p(n|Y, \tilde{\alpha}, \lambda) = \frac{ p(Y|n, \tilde{\alpha}) p(n|\lambda) }{ p(Y|\tilde{\alpha}, \lambda) }$$

$\tilde{\alpha}$$\lambda$

$$p(Y|\tilde{\alpha}, \lambda) = \sum_{n=1}^\infty p(Y|n, \tilde{\alpha}) p(n|\lambda)$$

$p(Y|n, \tilde{\alpha}) = 0$$n < m$

$$p(Y|\tilde{\alpha}, \lambda) = \frac{1}{(e^\lambda - 1)} \sum_{n=m}^\infty \frac{\Gamma(n\tilde{\alpha})\prod_{i=1}^n \Gamma(y_i + \tilde{\alpha})}{\Gamma(n\tilde{\alpha} + \sum_{i=1}^n y_i) \Gamma(\tilde{\alpha})^n} \cdot \frac{\lambda^n}{n!}$$

Menuju ke:

$$p(n|Y,\tilde{\alpha}, \lambda) = \frac{\Gamma(n\tilde{\alpha})\prod_{i=1}^n \Gamma(y_i + \tilde{\alpha})}{\Gamma(n\tilde{\alpha} + \sum_{i=1}^n y_i) \Gamma(\tilde{\alpha})^n} \cdot \frac{\lambda^n}{n!} \cdot \left(\sum_{j=m}^\infty \frac{\Gamma(j\tilde{\alpha})\prod_{i=1}^j \Gamma(y_i + \tilde{\alpha})}{\Gamma(j\tilde{\alpha} + \sum_{i=1}^j y_i) \Gamma(\tilde{\alpha})^j} \cdot \frac{\lambda^j}{j!}\right)^{-1}$$

$[m, \infty)$

Berikut ini contoh ceroboh dalam R:

logPosteriorN <- function(max, Y, lambda, alpha){
    m <- length(Y)
    sumy <- sum(Y)
    pp <- sapply(1:max, function(j){
        prior <- log(lambda)*j - log(exp(lambda)-1) - lgamma(j+1)
        posterior <- lgamma(alpha*j) + sum(lgamma(Y + alpha)) - j*lgamma(alpha) - lgamma(sumy + j*alpha)
        if( j > m ) { posterior <- posterior + (j-m)*lgamma(alpha) } 
        else if( j < m ) { posterior = -Inf }
        prior + posterior
        })
    evidence <- log(sum(exp(pp))) # there's no check that this converges
    pp - evidence
}

## with even representation of sides
Y <- c(10, 10, 10, 10)
post <- logPosteriorN(30, Y, 10, 1.2)
plot(1:30, exp(post), pch=19, type="b")

## with uneven representation of sides
Y <- c(1, 2, 1, 0, 0, 2, 1, 0, 1)
post <- logPosteriorN(30, Y, 10, 1.2)
plot(1:30, exp(post), pch=19, type="b")

$\tilde{\alpha}$$n$$\tilde{\alpha}$

Tentu saja, ini adalah salah satu pendekatan untuk estimasi. Anda akan segera menemukan orang lain (dari rasa Bayesian dan non-Bayesian) dengan sedikit pencarian.

Bagian II (Jawaban untuk berkomentar)

$Y = \{y_1, \dots, y_m, y_{m+1}, \dots, y_n \}$$\Omega = \{\omega_1, \dots, \omega_m, \omega_{m+1}, \dots, \omega_n\}$

$$\mathrm{Pr}(Y|\Omega, n) = \frac{\Gamma(\sum_{i=1}^n y_i + 1)}{\prod_{i=1}^n \Gamma(y_i + 1) } \prod_{i=1}^n \omega_i^{y_i}$$

$y \in \mathbb{N}$$y_1 \dots y_m > 0$$y_{m+1} \dots y_n = 0$$n$$n$

$$\mathrm{Pr}(n|\lambda) = \frac{\lambda^{n}}{(\exp\{\lambda\} - 1)n!},~n \in \mathbb{Z}^+$$

$\Omega$$\tilde{\alpha}$$n$

$$\mathrm{Pr}(\Omega|\tilde{\alpha}, n) = \frac{\Gamma(n\tilde{\alpha})}{\Gamma(\tilde{\alpha})^n} \prod_{i=1}^n \omega_i^{\tilde{\alpha}-1}$$

$$\mathrm{Pr}(Y|\tilde{\alpha}, n) = \int \mathrm{Pr}(Y|\Omega, n) \mathrm{Pr}(\Omega|\tilde{\alpha}, n) = \frac{\Gamma(n \tilde{\alpha})} {\Gamma(\sum_{i=1}^n y_i + n \tilde{\alpha}) \Gamma(\tilde{\alpha})^n} \prod_{i=1}^n \Gamma(y_i + \tilde{\alpha})$$

$n$$i \in \{1 \dots n\}$$j < i$$m \leq n$$Y$$n-m$$\mathcal{P}[Y]$

$n$

$$\mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, n) = \frac{n!}{(n-m)!} \mathrm{Pr}(Y|\tilde{\alpha}, n)$$

$n$

$$\mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, \lambda) = \sum_{j=m}^{\infty} \mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, n) \mathrm{Pr}(n|\lambda)$$

$$\mathrm{Pr}(n|\mathcal{P}[Y], \tilde{\alpha}, \lambda) = \frac{\mathrm{Pr}(\mathcal{P}[Y]|n, \tilde{\alpha}) \mathrm{Pr}(n|\lambda)}{\mathrm{Pr}(\mathcal{P}[Y]|\tilde{\alpha}, \lambda)}$$

Cukup colokkan dari definisi di atas. Sekali lagi, penyebutnya adalah seri tak terbatas yang akan menyatu dengan cepat: dalam model sederhana ini, MCMC tidak perlu memberikan perkiraan yang memadai.

Dengan memodifikasi kode R dari Bagian I:

logPosteriorN_2 <- function(max, Y, lambda, alpha){
    m <- length(Y)
    sumy <- sum(Y)
    pp <- sapply(1:max, function(j){
        prior <- log(lambda)*j - log(exp(lambda)-1) - lgamma(j+1)
        likelihood <- lchoose(j, m) + lgamma(m + 1) + lgamma(alpha*j) + sum(lgamma(Y + alpha)) - j*lgamma(alpha) - lgamma(sumy + j*alpha)
        if( j > m ) { likelihood <- likelihood + (j-m)*lgamma(alpha) } 
        else if( j < m ) { likelihood = -Inf }
        prior + likelihood
        })
    evidence <- log(sum(exp(pp))) # there's no check that this converges
    pp - evidence
}

Y_1 <- rep(10, 15)
pos_1 <- logPosteriorN_2(50, Y_1, 6, 1)
plot(1:50, exp(pos_1))