Description Usage Arguments Details Value Examples
AcceptProp
is a utility function to determine if a proposal should
be accepted in a Metropolis or Metropolis-Hastings step.
1 2 | AcceptProp(log.curr, log.prop, log.curr.to.prop = 0,
log.prop.to.curr = 0)
|
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')) |
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.
TRUE/FALSE
for whether the proposal should be accepted or
rejected, respectively
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 | # 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.