Application of Metropolis Hastings

Question

I am trying to implement the Metropolis Hastings algorithm for Bayesian analysis. In this case, the parameter of interest is the scale parameter for a Weibull distribution. The context is for reliability estimation; the scale parameter is analogous to the service life of an item.

I am working in R. I have data that has an "Age" property that is the age at which an item failed. There is no censored data. I fit the data to a Weibull distribution:

weib_fit <- fitdist(data_tbl$Age, "weibull")

I create a function to serve as the prior distribution based on some assumed 2-param Weibull distribution:

f <- function(x)
{
    shape = 3 # theoretical degradation
    scale = 30 # theoretical design life
    return(dweibull(x, shape, scale))
}

I create a liklihood function using the data parameters:

data_shape <- as.numeric(weib_fit$estimate["shape"])
data_scale <- as.numeric(weib_fit$estimate["scale"])

qx <- function(x)
{
    # replace the scale param with x
    rweibull(1, data_shape, x)
}

The rest is outlined below

step <- function(x, f, qx)
{
   xp <- qx(x)
   prob <- min(1, f(xp)/f(x))
   if (runif(1) < prob)
   {
        x <- xp
   }
   return(x)
}

run <- function(x, f, qx, nsteps)
{
    res <- matrix(NA, nsteps, length(x))
    for (i in seq_len(nsteps))
        res[i,] <- x <- step(x, f, qx)
    drop(res)
}

nsteps <- 100000
age_guess <- 20
res <- run(age_guess, f, qx, nsteps)

burn_in <- 5000
# get the final result
result <- res[c(burn_in:length(res))]
hist(result, main="MH Distribution", xlab="Service Life Param")

The problem: I don't think I'm understanding the implementation of the Metropolis-Hastings algorithm. I would assume the qx function above would be selecting for the Weibull scale distribution, but I get somewhat non-sensical results. If I modify it to rweibull(1, x, data_scale) it aligns much more with the expected frequentist results, but it doesn't make sense to me why I would be iterating by passing x to the rweibull() shape parameter when I am interested in the distribution of the scale parameter.

What am I missing in my understanding and how would I need to modify the code above to fix it?

Answer 1

There is a lot wrong with my question above that I came to realize after studying a bit more. I'll outline a different approach I've taken that works. The difference in this case is that I'm solving for a two parameter space instead of a single parameter space. The two parameters are the shape and scale parameter of a Weibull distribution for failures. This follows the blog entry:

https://theoreticalecology.wordpress.com/2010/09/17/metropolis-hastings-mcmc-in-r/

I define a likelihood function that accepts failure data and a two parameter array. The first index in the array aligns to the shape parameter and the second to the scale parameter. These are used to define the Weibull density function that serves as the likelihood function, that is used across all the individual failure data points.

likelihood <- function(data, params)
{
  single_likelihoods <- dweibull(data, shape=params[1], scale=params[2], log=TRUE)
  sumll <- sum(single_likelihoods)
  return(sumll)
}

I create a prior distribution. Here I use flat, uninformative priors. If I were using some other prior distribution, I may pass the function a 2-parameter params argument

prior <- function()
{
  shape_prior <- log(1)
  scale_prior <- log(1)
  
  return(shape_prior + scale_prior)
}

I create a posterior distribution which is the product of the likelihood and prior probabilities. (Since I'm use log values, this two probabilities are summed):

posterior <- function(data, params)
{
  return(likelihood(data, params) + prior())
}

I then need a function to create a proposed new shape/scale combination for each step. Here I assume the potential parameters are normally distributed.

proposalfunction <- function(param)
{
  proposed.params <- rnorm(length(param),mean = param, sd= c(.1, .5))
  return(proposed.params)
}

I then create a function to iterate over a specified number of steps:

MCMC_iter <- function(data, iterations){
  chain = array(dim = c(iterations+1,2))
  chain[1,] = c(2,30) # starting values for shape,scale
  
  for (i in 1:iterations){
    proposal = proposalfunction(chain[i,])
    
    probab = exp(posterior(data, proposal) - posterior(data, chain[i,]))
    if (runif(1) < probab){
      chain[i+1,] = proposal
    }else{
      chain[i+1,] = chain[i,]
    }
  }
  return(chain)
}

I can then create a markov chain using this function

n = 10000 # number of iterations
chain <- run_metropolis_MCMC(data_tbl$Age, n)

I define a burn-in:

burnin <- 5000
mcmc <- chain[burnin:n,]

To plot the posteriors for each post-burn-in:

hist(mcmc[,1]) #shape posterior
hist(mcmc[,2]) #scale posterior

To calculate the final shape/scale estimate:

  mcmc.shape <- mean(mcmc[,1])
  mcmc.scale <- mean(mcmc[,2])