TDD (version 0.1)

Metropolis: Metropolis-Hastings Markov Chain Monte Carlo

Description

Uses the Metropolis-Hastings Markov Chain Monte Carlo (MCMC) method to determine an optimal model to fit some data set.

Usage

Metropolis(loglikelihood, sigma, m1, niter, gen, logproposal, logprior = function(x) 0, burn = 0, save_int = 10, ...)

Value

List including the following elements:mMatrix where each row is the test model parameters of an iterationllog-likelihood of each iteration's modelaAcceptance ratio (moving window of length 100)bestList including best model found and its log-likelihood

Details

Metropolis prints progress information to the screen every 1000 iterations. These lines include the following: Number of iterations completed out of total loglikelihood of current model loglikelihood of proposed model loglikelihood of best model found so far Whether the proposed model this round is rejected or accepted Acceptance ratio over the last 100 iterations

References

Hastings, W.K. (1970). "Monte Carlo Sampling Methods Using Markov Chains and Their Applications". Biometrika 57 (1): 97-109.

Aster, R.C., Borchers, B., Thurber, C.H., Parameter Estimation and Inverse Problems, Elsevier Academic Press, 2012.

Examples

# try to find a non-negative vector m so that
# 1.  sqrt(m[1]^2 + m[2]^2) is close to 5
# 2.  sqrt(m[1] * m[2]) is close to 3.5
# 3.  2 * m[1] + m[2] is close to 10

# We are trying to match this data vector:
data = c(5, 3.5, 10)

# Define log-likelihood as -0.5 * sum of squared differences between
# modeled and true data
loglikelihood = function(model, data){
  d2 = c(sqrt(sum(model^2)), sqrt(abs(prod(model))), sum(model*c(2,1)))
  -0.5 * sum((d2 - data)^2)
}

# A proposed model is generated by randomly picking a model parameter
# and perturbing it by a random number distributed normally according to sigma
generate = function(x, sigma){
  w = ceiling(runif(1) * length(x))
  x[w] = x[w] + rnorm(1, 0, sigma[w])
  return(x)
}

# Proposal distribution is defined as multivariate normal, with mean
# zero and standard deviations sigma:
logproposal = function(x1, x2, sigma){
  -0.5 * sum(((x1) - (x2))^2/(sigma+1e-12)^2)
}

# logprior reflects prior knowledge that the answer is non-negative
logprior = function(m){
  if(all(m >= 0))
    0
  else
    -Inf
}


sigma = c(0.1, 0.1)
m1 = c(0, 0)
x = Metropolis(loglikelihood, sigma, m1, niter = 5000, gen = generate,
logproposal = logproposal, logprior = logprior, burn = 0, save_int = 1,
data = data)

# Notice the high acceptance ratios--this means that values in sigma are
# too low.  The MCMC is probably "optimally tuned" when sigma is set so 
# acceptance ratios vary between roughly 0.2 and 0.5.

# Plot models--par 1 on x, par 2 on y axis
# Note initial trajectory away from m1 (0, 0) to more likely
# region--this can be eliminated by setting 'burn' to a higher value
plot(x$m[,1], x$m[,2], pch = '.', col = rainbow(nrow(x$m)))

# Histograms/scatter plots showing posterior distributions.
# Note the strong covariance between these parameters!
par(mfrow = c(2, 2))
hist(x$m[,1])
plot(x$m[,2], x$m[,1], pch = '.')
plot(x$m[,1], x$m[,2], pch = '.')
hist(x$m[,2])