貝葉斯估計𝑁ñN二項分佈的
這個問題是這個問題的技術後續。
我無法理解和復制Raftery (1988) 中提出的模型：二項式推理參數：WinBUGS/OpenBUGS/JAGS 中的分層貝葉斯方法。它不僅僅是關於代碼的，所以它應該是這裡的主題。
背景
讓是一組來自未知的二項分佈的成功計數和. 此外，我假設遵循帶參數的泊松分佈（如論文中所述）。那麼，每個具有均值的泊松分佈. 我想指定先驗和.
假設我沒有任何好的先驗知識或者，我想為兩者分配非信息性先驗和. 說，我的先驗是和.
作者使用了不恰當的先驗但 WinBUGS 不接受不恰當的先驗。
例子
在論文（第 226 頁）中，提供了以下觀察到的 waterbucks 的成功計數：. 我要估計, 人口規模。
這是我嘗試在 WinBUGS 中解決示例的方法（在@Stéphane Laurent 的評論之後更新）：
model { # Likelihood for (i in 1:N) { x[i] ~ dbin(theta, n) } # Priors n ~ dpois(mu) lambda ~ dgamma(0.001, 0.001) theta ~ dunif(0, 1) mu <- lambda/theta } # Data list(x = c(53, 57, 66, 67, 72), N = 5) # Initial values list(n = 100, lambda = 100, theta = 0.5) list(n = 1000, lambda = 1000, theta = 0.8) list(n = 5000, lambda = 10, theta = 0.2)該模型在500,000個樣本和 20,000 個老化樣本後仍不能很好地收斂。這是 JAGS 運行的輸出：
Inference for Bugs model at "jags_model_binomial.txt", fit using jags, 5 chains, each with 5e+05 iterations (first 20000 discarded), n.thin = 5 n.sims = 480000 iterations saved mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff lambda 63.081 5.222 53.135 59.609 62.938 66.385 73.856 1.001 480000 mu 542.917 1040.975 91.322 147.231 231.805 462.539 3484.324 1.018 300 n 542.906 1040.762 95.000 147.000 231.000 462.000 3484.000 1.018 300 theta 0.292 0.185 0.018 0.136 0.272 0.428 0.668 1.018 300 deviance 34.907 1.554 33.633 33.859 34.354 35.376 39.213 1.001 43000問題
顯然，我遺漏了一些東西，但我看不出到底是什麼。我認為我對模型的表述在某個地方是錯誤的。所以我的問題是：
- 為什麼我的模型及其實現不起作用？
- Raftery (1988) 給出的模型如何正確地制定和實施？
謝謝你的幫助。
好吧，既然你的代碼可以工作，看起來這個答案有點太晚了。但是我已經寫了代碼，所以…
對於它的價值，這是相同的*模型適合
rstan. 在我的消費類筆記本電腦上估計需要 11 秒，為我們感興趣的參數實現更高的有效樣本量在更少的迭代中。raftery.model <- " data{ int I; int y[I]; } parameters{ real<lower=max(y)> N; simplex[2] theta; } transformed parameters{ } model{ vector[I] Pr_y; for(i in 1:I){ Pr_y[i] <- binomial_coefficient_log(N, y[i]) +multiply_log(y[i], theta[1]) +multiply_log((N-y[i]), theta[2]); } increment_log_prob(sum(Pr_y)); increment_log_prob(-log(N)); } " raft.data <- list(y=c(53,57,66,67,72), I=5) system.time(fit.test <- stan(model_code=raftery.model, data=raft.data,iter=10)) system.time(fit <- stan(fit=fit.test, data=raft.data,iter=10000,chains=5))請注意，我將其轉換
theta為 2-simplex。這只是為了數值穩定性。感興趣的數量是theta[1]; 顯然theta[2]是多餘的信息。*如你所見，後面的總結幾乎是相同的，並且促進實際數量似乎對我們的推論沒有實質性影響。
97.5% 的分位數比我的模型大 50%，但我認為這是因為 stan 的採樣器比簡單的隨機遊走更擅長探索後驗的全部範圍，因此它更容易進入尾部。不過，我可能是錯的。
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat N 1078.75 256.72 15159.79 94.44 148.28 230.61 461.63 4575.49 3487 1 theta[1] 0.29 0.00 0.19 0.01 0.14 0.27 0.42 0.67 2519 1 theta[2] 0.71 0.00 0.19 0.33 0.58 0.73 0.86 0.99 2519 1 lp__ -19.88 0.02 1.11 -22.89 -20.31 -19.54 -19.09 -18.82 3339 1取值從 stan 生成，我使用這些來繪製後驗預測值. 我們不應該對後驗預測的平均值感到驚訝非常接近樣本數據的平均值！
N.samples <- round(extract(fit, "N")[[1]]) theta.samples <- extract(fit, "theta")[[1]] y_pred <- rbinom(50000, size=N.samples, prob=theta.samples[,1]) mean(y_pred) Min. 1st Qu. Median Mean 3rd Qu. Max. 32.00 58.00 63.00 63.04 68.00 102.00為了檢查
rstan採樣器是否有問題，我計算了網格上的後驗。我們可以看到後部是香蕉形的；這種後驗對於歐幾里得度量 HMC 可能是有問題的。但是讓我們檢查一下數值結果。（香蕉形狀的嚴重性實際上在這裡被壓制了，因為是在對數刻度上。）如果你想一分鐘香蕉的形狀，你會意識到它一定在線上.
下面的代碼可以確認我們從 stan 得到的結果是有意義的。
theta <- seq(0+1e-10,1-1e-10, len=1e2) N <- round(seq(72, 5e5, len=1e5)); N[2]-N[1] grid <- expand.grid(N,theta) y <- c(53,57,66,67,72) raftery.prob <- function(x, z=y){ N <- x[1] theta <- x[2] exp(sum(dbinom(z, size=N, prob=theta, log=T)))/N } post <- matrix(apply(grid, 1, raftery.prob), nrow=length(N), ncol=length(theta),byrow=F) approx(y=N, x=cumsum(rowSums(post))/sum(rowSums(post)), xout=0.975) $x [1] 0.975 $y [1] 3236.665嗯。這不是我所期望的。97.5% 分位數的網格評估更接近 JAGS 結果而不是
rstan結果。同時，我不認為網格結果應該被視為福音，因為網格評估正在做一些相當粗略的簡化：一方面網格分辨率不是太精細，而另一方面，我們（錯誤地) 斷言網格中的總概率必須為 1，因為我們必須繪製邊界（和有限網格點）才能使問題可計算（我仍在等待無限 RAM）。事實上，我們的模型在. 但也許這裡有更微妙的東西在起作用。
