###########################
##Fit a spatial regression
###########################
set.seed(1)
n <- 50
x <- runif(n, 0, 100)
y <- runif(n, 0, 100)
D <- as.matrix(dist(cbind(x, y)))
phi <- 3/50
sigmasq <- 50
tausq <- 20
mu <- 150
s <- (sigmasq*exp(-phi*D))
w <- mvrnorm(1, rep(0, n), s)
Y <- mvrnorm(1, rep(mu, n) + w, tausq*diag(n))
X <- as.matrix(rep(1, length(Y)))
##Priors
##IG sigma^2 and tau^2
a.sig <- 2
b.sig <- 100
a.tau <- 2
b.tau <- 100
##Unif phi
a.phi <- 3/100
b.phi <- 3/1
##Functions used to transform phi to continuous support.
logit <- function(theta, a, b){log((theta-a)/(b-theta))}
logit.inv <- function(z, a, b){b-(b-a)/(1+exp(z))}
##Metrop. target
target <- function(theta){
mu.cand <- theta[1]
sigmasq.cand <- exp(theta[2])
tausq.cand <- exp(theta[3])
phi.cand <- logit.inv(theta[4], a.phi, b.phi)
Sigma <- sigmasq.cand*exp(-phi.cand*D)+tausq.cand*diag(n)
SigmaInv <- chol2inv(chol(Sigma))
logDetSigma <- determinant(Sigma, log=TRUE)$modulus[1]
out <- (
##Priors
-(a.sig+1)*log(sigmasq.cand) - b.sig/sigmasq.cand
-(a.tau+1)*log(tausq.cand) - b.tau/tausq.cand
##Jacobians
+log(sigmasq.cand) + log(tausq.cand)
+log(phi.cand - a.phi) + log(b.phi -phi.cand)
##Likelihood
-0.5*logDetSigma-0.5*(t(Y-X%*%mu.cand)%*%SigmaInv%*%(Y-X%*%mu.cand))
)
return(out)
}
##Run a couple chains
n.batch <- 500
batch.length <- 25
inits <- c(0, log(1), log(1), logit(3/10, a.phi, b.phi))
chain.1 <- adaptMetropGibbs(ltd=target, starting=inits,
batch=n.batch, batch.length=batch.length, report=100)
inits <- c(500, log(100), log(100), logit(3/90, a.phi, b.phi))
chain.2 <- adaptMetropGibbs(ltd=target, starting=inits,
batch=n.batch, batch.length=batch.length, report=100)
##Check out acceptance rate just for fun
plot(mcmc.list(mcmc(chain.1$acceptance), mcmc(chain.2$acceptance)))
##Back transform
chain.1$p.samples[,2] <- exp(chain.1$p.samples[,2])
chain.1$p.samples[,3] <- exp(chain.1$p.samples[,3])
chain.1$p.samples[,4] <- 3/logit.inv(chain.1$p.samples[,4], a.phi, b.phi)
chain.2$p.samples[,2] <- exp(chain.2$p.samples[,2])
chain.2$p.samples[,3] <- exp(chain.2$p.samples[,3])
chain.2$p.samples[,4] <- 3/logit.inv(chain.2$p.samples[,4], a.phi, b.phi)
par.names <- c("mu", "sigma.sq", "tau.sq", "effective range (-log(0.05)/phi)")
colnames(chain.1$p.samples) <- par.names
colnames(chain.2$p.samples) <- par.names
##Discard burn.in and plot and do some convergence diagnostics
chains <- mcmc.list(mcmc(chain.1$p.samples), mcmc(chain.2$p.samples))
plot(window(chains, start=as.integer(0.5*n.batch*batch.length)))
gelman.diag(chains)
##########################
##Example of fitting a
##a non-linear model
##########################
##Example of fitting a non-linear model
set.seed(1)
########################################################
##Simulate some data.
########################################################
a <- 0.1 #-Inf < a < Inf
b <- 0.1 #b > 0
c <- 0.2 #c > 0
tau.sq <- 0.1 #tau.sq > 0
fn <- function(a,b,c,x){
a+b*exp(x/c)
}
n <- 200
x <- seq(0,1,0.01)
y <- rnorm(length(x), fn(a,b,c,x), sqrt(tau.sq))
##check out your data
plot(x, y)
########################################################
##The log target density
########################################################
##Define the log target density used in the Metrop.
ltd <- function(theta){
##extract and transform as needed
a <- theta[1]
b <- exp(theta[2])
c <- exp(theta[3])
tau.sq <- exp(theta[4])
y.hat <- fn(a, b, c, x)
##likelihood
logl <- sum(dnorm(y, y.hat, sqrt(tau.sq), log=TRUE))
##priors IG on tau.sq and normal on a and transformed b, c, d
logl <- (logl
-(IG.a+1)*log(tau.sq)-IG.b/tau.sq
+sum(dnorm(theta[1:3], N.mu, N.v, log=TRUE))
)
##Jacobian adjustment for tau.sq
logl <- logl+log(tau.sq)
return(logl)
}
########################################################
##The rest
########################################################
##Priors
IG.a <- 2
IG.b <- 0.01
N.mu <- 0
N.v <- 10
theta.tuning <- c(0.01, 0.01, 0.005, 0.01)
##Run three chains with different starting values
n.batch <- 1000
batch.length <- 25
theta.starting <- c(0, log(0.01), log(0.6), log(0.01))
run.1 <- adaptMetropGibbs(ltd=ltd, starting=theta.starting, tuning=theta.tuning,
batch=n.batch, batch.length=batch.length, report=100)
theta.starting <- c(1.5, log(0.05), log(0.5), log(0.05))
run.2 <- adaptMetropGibbs(ltd=ltd, starting=theta.starting, tuning=theta.tuning,
batch=n.batch, batch.length=batch.length, report=100)
theta.starting <- c(-1.5, log(0.1), log(0.4), log(0.1))
run.3 <- adaptMetropGibbs(ltd=ltd, starting=theta.starting, tuning=theta.tuning,
batch=n.batch, batch.length=batch.length, report=100)
##Back transform
samples.1 <- cbind(run.1$p.samples[,1], exp(run.1$p.samples[,2:4]))
samples.2 <- cbind(run.2$p.samples[,1], exp(run.2$p.samples[,2:4]))
samples.3 <- cbind(run.3$p.samples[,1], exp(run.3$p.samples[,2:4]))
samples <- mcmc.list(mcmc(samples.1), mcmc(samples.2), mcmc(samples.3))
##Summary
plot(samples, density=FALSE)
gelman.plot(samples)
burn.in <- 5000
fn.pred <- function(theta,x){
a <- theta[1]
b <- theta[2]
c <- theta[3]
tau.sq <- theta[4]
rnorm(length(x), fn(a,b,c,x), sqrt(tau.sq))
}
post.curves <- t(apply(samples.1[burn.in:nrow(samples.1),], 1, fn.pred, x))
post.curves.quants <- summary(mcmc(post.curves))$quantiles
plot(x, y, pch=19, xlab="x", ylab="f(x)")
lines(x, post.curves.quants[,1], lty="dashed", col="blue")
lines(x, post.curves.quants[,3])
lines(x, post.curves.quants[,5], lty="dashed", col="blue")Run the code above in your browser using DataLab