overture (version 0.4-0)

AcceptProp: Determine if a Metropolis<U+2013>Hastings step should be accepted

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'))\)

Value

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

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.

Examples

# NOT RUN {
# 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)
# }