r - 试图从贝叶斯统计中复制数字而不流泪 : A sampling-resampling perspective, 但失败了

Question

我正在尝试不流泪地复制贝叶斯统计论文中的三个数字:采样-重采样视角，可在此处找到:http://hedibert.org/wp-content/uploads/2013/12/1992SmithGelfand.pdf我的目标是复制第 5 部分的结果。这是我的代码:

theta1<-runif(1000,0,1)
theta2<-runif(1000,0,1)
theta<-cbind(theta1,theta2)
theta<-as.data.frame(theta)

plot(theta1,theta2)

n1<-c(5,6,4)
n2<-c(5,4,6)
y<-c(7,5,6)

l<-rep(NA,nrow(theta))

for (i in 1:nrow(theta)){
  
  llh.1.store<-rep(NA,4)
  for (j in 2:5){
    llh.1.store[j-1]<-(factorial(n1[1])/(factorial(j)*factorial(n1[1]-j)))*(factorial(n2[1])/(factorial(y[1]-j)*factorial(n2[1]-y[1]+j)))*(theta[i,1]^j)*((1-theta[i,1])^(n1[1]-j))*(theta[i,2]^(y[1]-j))*((1-theta[i,2])^(n2[1]-y[1]+j))
  }
  llh1<-sum(llh.1.store)
  
  llh.2.store<-rep(NA,5)
  for (x in 1:5){
    llh.2.store[x]<-(factorial(n1[2])/(factorial(x)*factorial(n1[2]-x)))*(factorial(n2[2])/(factorial(y[2]-x)*factorial(n2[2]-y[2]+x)))*(theta[i,1]^x)*((1-theta[i,1])^(n1[2]-x))*(theta[i,2]^(y[2]-x))*((1-theta[i,2])^(n2[2]-y[2]+x))
  }
  llh2<-sum(llh.2.store)
  
  llh.3.store<-rep(NA,5)
  for (t in 0:4){
    llh.3.store[t+1]<-(factorial(n1[3])/(factorial(t)*factorial(n1[3]-t)))*(factorial(n2[3])/(factorial(y[3]-t)*factorial(n2[3]-y[3]+t)))*(theta[i,1]^t)*((1-theta[i,1])^(n1[3]-t))*(theta[i,2]^(y[3]-t))*((1-theta[i,2])^(n2[3]-y[3]+t))
  }
  llh3<-sum(llh.3.store)
  
  l[i]<-prod(llh1,llh2,llh3)
}

q<-l/sum(l)
post.theta<-sample_n(theta,1000,replace=TRUE,weight=q)

ggplot(post.theta) +
  aes(x = theta1, y = theta2) +
  geom_point(
    shape = "circle",
    size = 1.85,
    colour = "#440154"
  ) +
  labs(title = "Sample from Posterior") +
  ggthemes::theme_few()

但它不会生成与图 2 相同的图。任何人都可以告诉我我做错了什么吗？

Answer 1

下面是 Gelfand 和 Smith (1990, TAS) 模型的描述:

因此，可能性确实是由 j 的三个和的乘积构成的。

出于好奇，我为此后验编写了一个 Gibbs 采样器，将缺失的 X¹[i] 模拟为潜在变量，与 OP 中的上述结果相比，我没有发现结果(黄色)有任何差异(蓝色):

这是我的 R 代码:

T=1e4
    
theta=matrix(NA,nrow=T,ncol=2)
x1=rep(NA,3)
    
theta[1,]=runif(2)
for(t in 1:(T-1)){
   for(j in 1:3){
     a<-max(0,n2[j]-y[j]):min(y[j],n1[j])
     x1[j]=sample(a,1,
        prob=choose(n1[j],a)*choose(n2[j],y[j]-a)*
           theta[t,1]^a*(1-theta[t,1])^(n1[j]-a)*
           theta[t,2]^(y[j]-a)*(1-theta[t,2])^(n2[j]-y[j]+a)
           )
      }
      theta[t+1,1]=rbeta(1,sum(x1)+1,sum(n1)-sum(x1)+1)
      theta[t+1,2]=rbeta(1,sum(y)-sum(x1)+1,sum(n2)-sum(y)+sum(x1)+1)
}
    
plot(theta[sample(1:nrow(theta),1000),],
   xlab = expression(theta[1]), 
   ylab = expression(theta[2]),
   col = "#440154", pch=19, cex=.4,
   main = "Sample from Posterior")

这个结果也与从论文中得出的似然函数相符: