# NOT RUN {
# }
# NOT RUN {
## logistic regression with an improper uniform prior
## X and y are passed as args to MCMCmetrop1R
logitfun <- function(beta, y, X){
eta <- X %*% beta
p <- 1.0/(1.0+exp(-eta))
sum( y * log(p) + (1-y)*log(1-p) )
}
x1 <- rnorm(1000)
x2 <- rnorm(1000)
Xdata <- cbind(1,x1,x2)
p <- exp(.5 - x1 + x2)/(1+exp(.5 - x1 + x2))
yvector <- rbinom(1000, 1, p)
post.samp <- MCMCmetrop1R(logitfun, theta.init=c(0,0,0),
X=Xdata, y=yvector,
thin=1, mcmc=40000, burnin=500,
tune=c(1.5, 1.5, 1.5),
verbose=500, logfun=TRUE)
raftery.diag(post.samp)
plot(post.samp)
summary(post.samp)
## ##################################################
## negative binomial regression with an improper unform prior
## X and y are passed as args to MCMCmetrop1R
negbinfun <- function(theta, y, X){
k <- length(theta)
beta <- theta[1:(k-1)]
alpha <- exp(theta[k])
mu <- exp(X %*% beta)
log.like <- sum(
lgamma(y+alpha) - lfactorial(y) - lgamma(alpha) +
alpha * log(alpha/(alpha+mu)) +
y * log(mu/(alpha+mu))
)
}
n <- 1000
x1 <- rnorm(n)
x2 <- rnorm(n)
XX <- cbind(1,x1,x2)
mu <- exp(1.5+x1+2*x2)*rgamma(n,1)
yy <- rpois(n, mu)
post.samp <- MCMCmetrop1R(negbinfun, theta.init=c(0,0,0,0), y=yy, X=XX,
thin=1, mcmc=35000, burnin=1000,
tune=1.5, verbose=500, logfun=TRUE,
seed=list(NA,1))
raftery.diag(post.samp)
plot(post.samp)
summary(post.samp)
## ##################################################
## sample from a univariate normal distribution with
## mean 5 and standard deviation 0.1
##
## (MCMC obviously not necessary here and this should
## really be done with the logdensity for better
## numerical accuracy-- this is just an illustration of how
## MCMCmetrop1R works with a density rather than logdensity)
post.samp <- MCMCmetrop1R(dnorm, theta.init=5.3, mean=5, sd=0.1,
thin=1, mcmc=50000, burnin=500,
tune=2.0, verbose=5000, logfun=FALSE)
summary(post.samp)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab