mtrp: Metropolis-Hastings Algorithm Using Normal Distribution

In hjy78/mcmcn: Markov Chain Monte Carlo Sampler

Description

MCMC method for continuous random vector with normal distribution as proposal distribution using a random walk Metropolis algorithm.

Usage

mtrp(f, n, init, stepsize = 1, burn = 1000)

Arguments

f	Density function from which one wants to sample.
n	The numbers of samples one wants to obtain.
init	The initial value vector, which indicates the dimensions.
stepsize	A vector with stepsize[i] being the standard deviation of proposal for updating 'i'th variable.
burn	Times of iterations one wants to omit before recording.

Value

A "mcmcn" object 'list("chain" = chain, "reject" = k/iters, "acpt" = acpt)' with chain storing samples by row, reject being the rejection rate, acpt being whether to be accepted each iter.

Examples

# Generating Multivariate Normal Distribution Samples------------------------

# multivariate normal distribution pdf(mu, A)
# mu <- c(1, 3)
# A <- matrix(c(1, 0.1, 0.1, 1), nrow = 2)
f <- pdff("norm", c(1, 3), matrix(c(1, 0.1, 0.1, 1), nrow = 2))

# generating random variates using function `mtrp`
x.norm <- mtrp(f, 10000, c(3, 3), burn = 0)

# exploring the results
summary(x.norm)
plot(x.norm)
qqnorm(x.norm$chain[, 1], main = "QQ plot, 1st variable")

# Exploring the Results of Different Initial Values--------------------------

# exploring the results of different initial values
set.seed(1234)
par(mfrow = c(2, 1))
x.norm1 <- mtrp(f, 10000, c(-4, 3), burn = 0)
x.norm2 <- mtrp(f, 10000, c(1, 3), burn = 0)
x.norm3 <- mtrp(f, 10000, c(6, 3), burn = 0)
plot(1:500, x.norm1$chain[1:500, 1], type = "l", col = "red",
     ylab = "1st variable", ylim = c(-4, 6), main = "first 500 iters")
lines(1:500, x.norm2$chain[1:500, 1], col = "green")
lines(1:500, x.norm3$chain[1:500, 1], col = "blue")
points(rep(1, 3), c(-4, 1, 6), col = c("red", "green", "blue"), pch = 20)
plot(9501:10000, x.norm1$chain[9501:10000, 1], type = "l", col = "red",
     ylab = "1st variable", ylim = c(-4, 6), main = "last 500 iters")
lines(9501:10000, x.norm2$chain[9501:10000, 1], col = "green")
lines(9501:10000, x.norm3$chain[9501:10000, 1], col = "blue")
par(mfrow = c(1, 1))

# Exploring the Results of Different Stepsizes-------------------------------

# exploring the results of different stepsizes
set.seed(1234)
par(mfrow = c(2, 1))
for (stepsize in exp(seq(-2, 1, length.out = 4))) {
  x.norm <- mtrp(f, 1000, c(5, 3), stepsize = stepsize, burn = 0)
  plot(
    x.norm$chain[, 1],
    type = "l", ylab = "1st variable",
    main = paste("stepsize", round(stepsize,2), "rejection rate", round(x.norm$reject,2))
  )
}
par(mfrow = c(1, 1))

hjy78/mcmcn documentation built on Jan. 1, 2020, 1:03 p.m.