好吧，既然你的代码可以工作，看起来这个答案有点太晚了。但是我已经写了代码，所以...

对于它的价值，这是相同的*模型适合rstan. 在我的消费类笔记本电脑上估计需要 11 秒，为我们感兴趣的参数实现了更高的有效样本量$(N, \theta)$在更少的迭代中。

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]是多余的信息。

*如你所见，后面的总结几乎是相同的，并且促进$N$实际数量似乎对我们的推论没有实质性影响。

97.5% 的分位数$N$比我的模型大 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

取值$N, \theta$从 stan 生成，我使用这些来绘制后验预测值$\tilde{y}$. 我们不应该对后验预测的平均值感到惊讶$\tilde{y}$非常接近样本数据的平均值！

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 可能是有问题的。但是让我们检查一下数值结果。（香蕉形状的严重性实际上在这里被压制了，因为$N$是在对数刻度上。）如果你想一分钟香蕉的形状，你会意识到它一定在线上$\bar{y}=\theta N$.

下面的代码可以确认我们从 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）。事实上，我们的模型在$(0,1)\times\left\{N|N\in\mathbb{Z}\land N\ge72)\right\}$. 但也许这里有更微妙的东西在起作用。

再次感谢 @StéphaneLaurent 和 @user777 在评论中提供的宝贵意见。在对之前的进行一些调整之后$\lambda$我现在可以复制 Raftery (1988) 论文的结果。

这是我使用 JAGS 和 R 的分析脚本和结果：

#===============================================================================================================
# Load packages
#===============================================================================================================

sapply(c("ggplot2"
         , "rjags"
         , "R2jags"
         , "hdrcde"
         , "runjags"
         , "mcmcplots"
         , "KernSmooth"), library, character.only = TRUE)

#===============================================================================================================
# Model file
#===============================================================================================================

cat("
    model {

    # Likelihood    
    for (i in 1:N) {
      x[i] ~ dbin(theta, n)
    }

    # Prior       
    n ~ dpois(mu)
    lambda ~ dgamma(0.005, 0.005)
#     lambda ~ dunif(0, 1000)
    mu <- lambda/theta
    theta ~ dunif(0, 1)    
}    
", file="jags_model_binomial.txt")


#===============================================================================================================
# Data
#===============================================================================================================

data.list <- list(x = c(53, 57, 66, 67, 72, NA), N = 6) # Waterbuck example from Raftery (1988)

#===============================================================================================================
# Inits
#===============================================================================================================

jags.inits <- function() { 
  list(
    n = sample(max(data.list$x, na.rm = TRUE):1000, size = 1) 
    , theta = runif(1, 0, 1)
    , lambda = runif(1, 1, 10)
#     , cauchy  = runif(1, 1, 1000)
    #     , mu = runif(1, 0, 5)
  )
}

#===============================================================================================================
# Run the chains
#===============================================================================================================

# Parameters to store

params <- c("n"
            , "theta"
            , "lambda"
            , "mu"
            , paste("x[", which(is.na(data.list[["x"]])), "]", sep = "")
)

# MCMC settings

niter <- 500000 # number of iterations
nburn <- 20000  # number of iterations to discard (the burn-in-period)
nchains <- 5    # number of chains

# Run JAGS

out <- jags(
  data                 = data.list
  , parameters.to.save = params
  , model.file         = "jags_model_binomial.txt"
  , n.chains           = nchains
  , n.iter             = niter
  , n.burnin           = nburn
  , n.thin             = 50
  , inits              = jags.inits
  , progress.bar       = "text")

在我的台式电脑上计算大约需要 98 秒。

#===============================================================================================================
# Inspect results
#===============================================================================================================

print(out
      , digits = 2
      , intervals = c(0.025, 0.1, 0.25, 0.5, 0.75, 0.9,  0.975))

结果是：

Inference for Bugs model at "jags_model_binomial.txt", fit using jags,
 5 chains, each with 5e+05 iterations (first 20000 discarded), n.thin = 50
 n.sims = 48000 iterations saved
         mu.vect sd.vect  2.5%    10%    25%    50%    75%     90%   97.5% Rhat n.eff
lambda     62.90    5.18 53.09  56.47  59.45  62.74  66.19   69.49   73.49    1 48000
mu        521.28  968.41 92.31 113.02 148.00 232.87 467.10 1058.17 3014.82    1  1600
n         521.73  968.54 95.00 114.00 148.00 233.00 467.00 1060.10 3028.00    1  1600
theta       0.29    0.18  0.02   0.06   0.13   0.27   0.42    0.55    0.66    1  1600
x[6]       63.03    7.33 49.00  54.00  58.00  63.00  68.00   72.00   78.00    1 36000
deviance   34.88    1.53 33.63  33.70  33.85  34.34  35.34   36.81   39.07    1 48000

的后验均值$N$是$522$并且后中位数是$233$. 我计算了后验模式$N$在对数尺度上并对估计值进行反向转换：

jagsfit.mcmc <- as.mcmc(out)
jagsfit.mcmc <- combine.mcmc(jagsfit.mcmc)

hpd.80 <- hdr.den(log(as.vector(jagsfit.mcmc[, "n"])), prob = c(80), den = bkde(log(as.vector(jagsfit.mcmc[, "n"])), gridsize = 10000))

exp(hpd.80$mode)

[1] 149.8161

以及 80%-HPD$N$是：

(hpd.ints <- HPDinterval(jagsfit.mcmc, prob = c(0.8)))

               lower      upper
deviance 33.61011007  35.677810
lambda   56.08842502  69.089507
mu       72.42307587 580.027182
n        78.00000000 578.000000
theta     0.01026193   0.465714
x[6]     53.00000000  71.000000

后验模式为$N$是$150$并且 80%​​-HPD 是$(78; 578)$这非常接近论文中给出的限制，即$(80; 598)$.