mtrp_expu: Metropolis-Hastings Algorithm Using Exponential-Uniform...

In hjy78/mcmcn: Markov Chain Monte Carlo Sampler

Description

MCMC method for continuous random vector with exponential-uniform mixture distribution as proposal distribution using Metropolis-Hastings algorithm.

Usage

mtrp_expu(f, n, init, a = 0, rate = 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.
a	The left support boundary vector. if a numeric number is provided, it will be broadcast to a vector with the same length as init.
rate	A vector with rate[i] being the rate of proposal for updating 'i'th variable. if a numeric number is provided, it will be broadcast to a vector with the same length as init.
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

# the Exponential-Uniform Mixture--------------------------------------------

# suppose X ~ exponential-uniform mixture with parameters(a = 0, b = 1, rate = 1)
# the plot below shows the different pdfs of X and Exp(1)
curve(0.5 * exp(1-x), from = 1, xlim = c(0, 3), ylim = c(0, 1), ylab = "y")
lines(c(0, 1), c(0.5, 0.5))
curve(exp(-x), col = "red", add = TRUE)
legend("topright", c("Exp(1)", "X"), lty = c(1, 1), col = c("red", "black"))

# the Exponential-Gamma mixture----------------------------------------------
# (Suppose that the rate parameter Λ has Gamma ( r, beta ) distribution and X has
# Exp(Λ) distribution. That is, X|Λ = λ ∼ f_X (x|λ) = λe^{−λy}.)

# pdf f(x) = r * beta^r / (x + beta)^{r+1}
pareto2d <- function(r, beta) {
  function(x) ifelse(x[1] >= 0 && x[2] >=0, x[2]^r * exp(-x[2] * (x[1]+beta)), 0)
}
r <- 1
beta <- 2
f <- pareto2d(r, beta)

# generating random variates using function `mtrp_expu`
x.pareto <- mtrp_expu(f, 10000, 1, rate = 0.5, burn = 0)

# exploring the results
summary(x.pareto0)
plot(x.pareto0)
hist(x.pareto0$chain[, 1], freq = FALSE, main = "Histogram of Samples", xlab = "X")
curve(r * beta^r * (x + beta)^(-(r+1)), from = 0, col = "red", add = TRUE)

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