AcceptProp: Determine if a Metropolis–Hastings step should be accepted

In overture: Tools for Writing MCMC

Description

AcceptProp is a utility function to determine if a proposal should be accepted in a Metropolis or Metropolis-Hastings step.

Usage

AcceptProp(log.curr, log.prop, log.curr.to.prop = 0,
  log.prop.to.curr = 0)

Arguments

log.curr	log density of the target at the current value, log(P(x))
log.prop	log density of the target at the proposed value, log(P(x'))
log.curr.to.prop	log of transition distribution from current value to proposed value, log(g(x'|x))
log.prop.to.curr	log of transition distribution from proposed value to current value, log(g(x|x'))

Details

The function uses the Metropolis choice for a Metropolis/Metropolis-Hastings sampler, which accepts a proposed value x' with probability

A(x', x) = min(1, P(x')/P(x) g(x|x')/g(x'|x))

where P(x) is the target distribution and g(x'|x) is the proposal distribution.

Value

TRUE/FALSE for whether the proposal should be accepted or rejected, respectively

Examples

# Sample from triangular distribution P(x) = -2x + 2 ----------------------
# Target distribution
LogP <- function(x) {
    log(-2*x + 2)
}

# Generate proposals using Beta(1/2, 1/2)
shape1 <- 1/2
shape2 <- 1/2

RProp <- function() { # Draw proposal
    rbeta(1, shape1, shape2)
}

DLogProp <- function(x) { # Log density of proposal distribution
    dbeta(x, shape1, shape2, log=TRUE)
}

SampleX <- function(x) { # Draw once from the target distribution
    x.prop <- RProp()
    if(AcceptProp(LogP(x), LogP(x.prop), DLogProp(x.prop), DLogProp(x))) {
        x <- x.prop
    }

    return(x)
}

# Draw from the target distribution
n.samples <- 10000
samples <- vector(length=n.samples)
x <- 0.5
Mcmc <- InitMcmc(n.samples)
samples <- Mcmc({
    x <- SampleX(x)
})

# Plot the results
hist(samples$x, freq=FALSE, ylim=c(0, 2.5), xlim=c(0, 1), xlab="x")
grid <- seq(0, 1, length.out=500)
lines(grid, exp(LogP(grid)), col="blue")
legend("topright", legend="True density", lty=1, col="blue", cex=0.75)

overture documentation built on Aug. 11, 2019, 1:04 a.m.