exponentialdecay.R

albert — Oct 14, 2013, 10:43 AM

# Nonlinear regression problem
# See http://www.dme.ufrj.br/mcmc/exponentialdecaymodel-R2WinBUGS.html

# Illustrate fitting an exponential decay regression model
# using both LearnBayes and JAGS

######################################################
# Using LearnBayes
######################################################

# Load in the LearnBayes package

library(LearnBayes)

# Write a function defining the log posterior
# parameters are (beta1, beta1, log(sigma))

expdecay <- function(theta, data){
  beta1 <- theta[1]
  beta2 <- theta[2]
  sigma <- exp(theta[3])
  y <- data$y
  x <- data$x
  log.like <- sum(dnorm(y, 
                        beta1 * (1 - exp(-beta2 * x)), 
                        sigma, log=TRUE))
  return(log.like)
}

# read in data as a list

b.data <- list(x=c(1,2,3,4,5,7),
               y=c(8.3,10.3,19,16,15.6,19.8))

# find posterior mode and approx var/cov matrix using laplace function

l.fit <- laplace(expdecay, c(2, 2, 1), b.data)

# construct a random walk Metropolis chain
# 10,000 iterations using a proposal var-cov from laplace
# and scale = 2

rw.fit <- rwmetrop(expdecay,
                   list(var=l.fit$var, scale=2),
                   c(20, 0.4, 1),
                   100000, b.data)

# check the simulated draws

library(coda)

Loading required package: lattice

plot(mcmc(rw.fit$par))

# construct a smoothed scatterplot of the simulated draws of beta1 and beta2

with(rw.fit, smoothScatter(par[, 1:2]))

KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009

#############################################
# USING JAGS
#############################################

# Load in the rjags package

library(rjags)

Linked to JAGS 3.3.0
Loaded modules: basemod,bugs

# Write model description (note subtle change in model)

modelString = "
model{
for (i in 1:N){
  mu[i] <- beta1*(1-exp(-beta2*x[i]))
  y[i]~dnorm(mu[i],tau2)
}
beta1~dunif(-20,50)
beta2~dunif(-2,6)
tau2 ~ dgamma(0.001,0.001)
}
"
writeLines(modelString, con="expdecay.bug")

# Define a list containing the data and prior parameters

library(LearnBayes)

b.data <- list(x=c(1,2,3,4,5,7),
               y=c(8.3,10.3,19,16,15.6,19.8),
               N=6)

# Define initial values for chain 

b.inits <- list(beta1=18, beta2=1.78, tau2=1)

# Run the function jags.model to initialize the MCMC sampler

jags <- jags.model('expdecay.bug',
                   data = b.data,
                   inits = b.inits,
                   n.chains = 1,
                   n.adapt = 100)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 49

Initializing model

# Run the MCMC 1000 iterations as a burn-in

update(jags, 1000)

# Collect the simulated draws of the proportion for 10,000 iterations

posterior <- coda.samples(jags, c("beta1", "beta2", "tau2"),
                          n.iter=10000)

# Look at simulated draws

plot(posterior)

# Look at summaries

summary(posterior)

Iterations = 1101:11100
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

        Mean    SD Naive SE Time-series SE
beta1 18.919 4.372  0.04372        0.19499
beta2  1.107 1.201  0.01201        0.04230
tau2   0.135 0.104  0.00104        0.00215

2. Quantiles for each variable:

         2.5%     25%    50%    75%  97.5%
beta1 12.6445 16.4789 18.267 20.342 30.936
beta2  0.1894  0.4642  0.652  1.071  5.065
tau2   0.0158  0.0591  0.108  0.184  0.399

# Collect simulated draws of (beta1, beta2, and sigma)

BETA1 <- posterior[[1]][, "beta1"]
BETA2 <- posterior[[1]][, "beta2"]
SIGMA <- sqrt(1/posterior[[1]][, "tau2"])

# Plot data and overlay fitted model (notice I use medians to summarize)

with(b.data, plot(x, y))
curve(median(BETA1) - median(BETA1) * exp(-median(BETA2) * x),
                                      add=TRUE)