最近，看了一篇关于贝叶斯的文章：Bayesian Statistics without tears:A Sampling-Resampling Perspective.（http://faculty.chicagobooth.edu/hedibert.lopes/teaching/ccis2010/1992SmithGelfand.pdf）基本上看懂了“抽样-在抽样”技术。但是，文中第五部分例子结果和我自己做的不一样。不知哪位大侠看过本文，而且做出了上面的那个例子？求指教！
我的代码是这样的：
<br />
n <- matrix(c(5, 6, 4,<br />
5, 4, 6), nrow = 3)<br />
y <- c(7, 5, 6)</p>
<p>fun <- function(i, theta1, theta2) {<br />
fun2 <- function(ji, i, theta1, theta2) {</p>
<p> ncol(combn(n[i, 1], ji)) * ncol(combn(n[i, 2], y[i] - ji))*<br />
theta1^ji * (1-theta1)^(n[i, 1]-ji) *<br />
theta2^(y[i] - ji) * (1-theta2)^(n[i, 2]-y[i]+ji)</p>
<p> }<br />
lowbond <- max(c(0, y[i]-n[i, 2]))<br />
upperbond <- min(c(n[i, 1], y[i]))<br />
sum(sapply(lowbond:upperbond, fun2, i, theta1, theta2))</p>
<p>}</p>
<p>likelihood <- function(theta1, theta2) { #定义似然函数<br />
prod(sapply(1:3, fun, theta1, theta2))<br />
}</p>
<p>set.seed(123)<br />
theta1 <- runif(500)<br />
theta2 <- runif(500)<br />
plot(theta1, theta2, cex = 0.2)<br />
weight <- NULL<br />
for(i in 1:500) {<br />
weight <- c(weight, likelihood(theta1[i], theta2[i]))<br />
}<br />
windows()<br />
inde <- sample(1:500, prob=weight, replace = T) #重抽样<br />
plot(jitter(theta1[inde]), jitter(theta2[inde]), cex = 0.2) #稍微抖动一下<br />
</p>
