骰子有幾個面？JAGS 中的貝葉斯推理
問題
我想對一個類似於死邊數未知的系統做一些推斷。骰子滾動了幾次，之後我想推斷一個參數的概率分佈，該參數對應於骰子的面數，θ。
直覺
如果在 40 次滾動後您觀察到 10 個紅色、10 個藍色、10 個綠色和 10 個黃色，則 θ 似乎應該在 4 處達到峰值，並且滾動每一側的偏差是以 1/4 為中心的分佈。
θ 有一個微不足道的下限，即在數據中觀察到的不同邊的數量。
上限仍然未知。可能存在第五面，其偏差可能較低。您觀察到的缺少第五類的數據越多，θ = 4 的後驗概率就越高。
方法
我已經將 JAGS 用於類似的問題（通過 R 和 rjags），這在此處似乎是合適的。
關於數據，可以說
obs <- c(10, 10, 10, 10)對應於上面示例中的觀察結果。我認為觀察結果應該用多項分佈來建模
obs ~ dmulti(p, n)，其中p ~ ddirch(alpha)和n <- length(obs)。θ 與 隱含的類別數量相關聯
alpha，那麼我如何建模alpha以包含不同可能的類別數量？備擇方案？
我對貝葉斯分析還很陌生，所以可能完全是錯誤的樹，是否有替代模型可以為這個問題提供不同的見解？
非常感謝！大衛
這是一個有趣的問題，稱為“物種抽樣”，多年來受到了很多關注，並且包含許多其他估計問題（例如標記重新捕獲）。可以這麼說，在這種情況下，JAGS 對您沒有幫助——JAGS 無法處理迭代中具有可變維度的馬爾可夫鏈。必須求助於為此類問題設計的 MCMC 方案，例如可逆跳躍 MCMC。
這是一種適合您所描述的特定模型的方法，這是我在 Jeff Miller 的工作中第一次遇到的 ( arxived )。
第一部分（原始問題）
我將做出的一個假設是，對給定類別的觀察意味著存在較低等級的類別。也就是說，觀察第 9 面的擲骰子意味著存在 1-8 面。不一定要這樣——類別可以是任意的——但我會在我的例子中假設。這意味著與其他物種估計問題相比，0 值是可觀察的。
假設我們有一個多項式樣本，
在哪裡是觀察到的最大類別，是（未知）類別數，並且所有等於 0. 參數是有限的，我們需要先驗。任何離散的，適當的先驗，支持將工作; 以零截斷泊松為例：
多項式概率的一個方便的先驗是 Dirichlet，
並且簡單地假設.
為了使問題更易於處理，我們將權重邊緣化：
在這種情況下，它導致了經過充分研究的Dirichlet 多項式分佈。然後目標是估計條件後驗，
我明確假設的地方和是固定的超參數。很容易看出：
在哪裡在哪裡. 這個無限級數應該很快收斂（只要先驗的尾部不太重），因此很容易近似。對於截斷的泊松，它具有以下形式：
導致：
哪個支持. 在這種情況下不需要 MCMC，因為貝葉斯規則分母中的無限級數可以在不費力的情況下進行近似。
這是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")您的直覺是正確的：跨類別的稀疏抽樣會導致類別總數的不確定性更大。如果你想治療作為未知參數，您需要使用 MCMC 和和.
當然，這是一種估計方法。稍加搜索，您將很容易找到其他（貝葉斯和非貝葉斯風味）。
第二部分（評論回答）
是具有相應概率的部分觀測多項式向量:
在哪裡,和但除此之外，指數是任意的。和以前一樣，問題是推斷類別的真實數量, 我們從先驗開始例如零截斷泊松：
和以前一樣，我們處理多項式概率作為具有對稱超參數的狄利克雷分佈，即對於給定的,
對概率向量進行積分（邊緣化）得到多項 Dirichlet：
這是我們與上面第一部分中的模型不同的地方。在第一部分，有一個隱含的類別排序：例如，在一個-side die，類別（邊）有一個隱含的順序和任何類別的觀察意味著存在較小的類別. 在第二部分中，我們有一個部分觀察到的多項式隨機向量，它沒有隱式排序。換句話說，數據表示數據點的無序劃分為觀察到的類別。我將表示由以下結果產生的無序分區增強未觀察到的類別，如.
以真實類別數量為條件的無序分區的概率, 可以通過考慮導致相同分區的類別排列的數量來找到：
這可以集成在給：
使用貝葉斯規則檢索後驗：
只需從上面的定義中插入即可。同樣，分母是一個無限級數，會很快收斂：在這個簡單的模型中，MCMC 不需要給出適當的近似值。
通過修改第一部分的 R 代碼：
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))