x <- out$x
y <- out$y
lupost <- function(beta, x, y) {
eta <- as.numeric(x %*% beta)
logp <- ifelse(eta < 0, eta - log1p(exp(eta)), -log1p(exp(-eta)))
logq <- ifelse(eta < 0, -log1p(exp(eta)), -eta - log1p(exp(-eta)))
logl <- sum(logp[y == 1]) + sum(logq[y == 0])
return(logl - sum(beta^2)/8)
}
在完成上面数据以及函数的定义(它们都是mcmc算法的输入参数)之后，我们就可以利用mcmc包中的metrop()来产生随机数据，模拟模型的后验分布。也就是说，我们将要从后验分布中进行取样。
set.seed(42)
beta.init <- as.numeric(coefficients(out))
out <- metrop(lupost, beta.init, 1000, x = x, y = y)
names(out)
扫码加好友，拉您进群
全部版块
我的主页
收藏
