Metropolis-Hasting

Implementation of the Metropolis-Hasting sampler for Bayesian inference.

MH(logposterior, proposal, init_est, d)

Construct a Sampler object for Metropolis-Hasting sampling.

Arguments

  • logposterior : the logposterior of the parameter of interest.
  • proposal : the proposal distribution for random steps of the MCMC.
  • init_est : the initial/starting value for the markov chain.
  • d : the dimension of the posterior distribution.

Value

Returns a MH type object.

Example

In order to illustrate the modeling, the data is simulated from a simple linear regression expectation function. That is the model is given by

y_i = w_0 + w_1 x_i + e_i,   e_i ~ N(0, 1 / a)

To do so, let B = [w_0, w_1]' = [.2, -.9]', a = 1 / 5. Generate 200 hypothetical data:

library(gridExtra)
library(lattice)
library(StochMCMC)

set.seed(123)

# Define data parameters
w0 <- -.3; w1 <- -.5; stdev <- 5.; a <- 1 / stdev

# Generate Hypothetical Data
n <- 200;
x <- runif(n, -1, 1);
A <- cbind(1, x);
B <- rbind(w0, w1);
f <- A %*% B;
y <- f + rnorm(n, 0, a);

my_df = data.frame(Independent = round(x, 4), Dependent = round(y, 4));

To view the head of the data, run the following:

head(my_df)
#   Independent Dependent
# 1     -0.4248   -0.2297
# 2      0.5766   -0.5369
# 3     -0.1820   -0.2583
# 4      0.7660   -0.7525
# 5      0.8809   -0.9308
# 6     -0.9089    0.1454

Next is to plot this data which can be done as follows:

xyplot(Dependent ~ Independent, data = my_df, type = c("p", "g"), col = "black")
alternate text

In order to proceed with the Bayesian inference, the parameters of the model is considered to be random modeled by a standard Gaussian distribution. That is, B ~ N(0, I), where 0 is the zero vector. The likelihood of the data is given by,

L(w|[x, y], b) = ∏_{i=1}^n N([x_i, y_i]|w, b)

Thus the posterior is given by,

P(w|[x, y]) ∝ P(w)L(w|[x, y], b)

To start programming, define the probabilities

# The log prior function is given by the following codes:
logprior <- function(theta, mu = zero_vec, s = eye_mat) {
    w0_prior <- dnorm(theta[1], mu[1], s[1, 1], log = TRUE)
    w1_prior <- dnorm(theta[2], mu[2], s[2, 2], log = TRUE)
    w_prior <- c(w0_prior, w1_prior)

    w_prior %>% sum %>% return
}

# The log likelihood function is given by the following codes:
loglike <- function(theta, alpha = a) {
    yhat <- theta[1] + theta[2] * x

    likhood <- numeric()
    for (i in 1:length(yhat)) {
        likhood[i] <- dnorm(y[i], yhat[i], alpha, log = TRUE)
    }

    likhood %>% sum %>% return
}

# The log posterior function is given by the following codes:
logpost <- function(theta) {
    loglike(theta, alpha = a) + logprior(theta, mu = zero_vec, s = eye_mat)
}

To start the estimation, define the necessary parameters for the Metropolis-Hasting algorithm

# Hyperparameters
zero_vec <- c(0, 0)
eye_mat <- diag(2)

Run the MCMC:

set.seed(123);
mh_object <- MH(logpost, init_est = c(0, 0))
chain1 <- mcmc(mh_object, r = 10000)

Extract the estimate

burn_in <- 100;
thinning <- 10;

# Expetation of the Posterior
est1 <- colMeans(chain1[seq((burn_in + 1), nrow(chain1), by = thinning), ])
est1
# [1] -0.2984246 -0.4964463