Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

In answering yet another question on X validated about the numerical approximation of the marginal likelihood, I suggested using an harmonic mean estimate as a simple but worthless solution based on an MCMC posterior sample. This was on a toy example with a uniform prior on (0,π) and a “likelihood” equal to sin(θ) [really a toy problem!]. Simulating an MCMC chain by a random walk Metropolis-Hastings algorithm is straightforward, as is returning the harmonic mean of the sin(θ)’s.

f <- function(x){
if ((0<x)&(x<pi)){
return(sin(x))}else{
return(0)}}

n = 2000 #number of iterations
sigma = 0.5
x = runif(1,0,pi) #initial x value
chain = fx = f(x)
#generates an array of random x values from norm distribution
rands = rnorm(n,0, sigma)
#Metropolis - Hastings algorithm
for (i in 2:n){
can = x + rands[i]  #candidate for jump
fcan=f(can)
aprob = fcan/fx #acceptance probability
if (runif(1) < aprob){
x = can
fx = fcan}
chain=c(chain,fx)}
I = pi*length(chain)/sum(1/chain) #integral harmonic approximation

However, the outcome looks remarkably stable and close to the expected value 2/π, despite 1/sin(θ) having an infinite integral on (0,π). Meaning that the average of the 1/sin(θ)’s has no variance. Hence I wonder why this specific example does not lead to an unreliable output… But re-running the chain with a smaller scale σ starts producing values of sin(θ) regularly closer to zero, which leads to an estimate of I both farther away from 2 and much more variable. No miracle, in the end!

Filed under: Books, Kids, Mountains, pictures, R, Running, Statistics, Travel Tagged: Gaussian random walk, harmonic mean estimator, Metropolis-Hastings algorithm, Monte Carlo Statistical Methods, numerical integration, simulation