# puzzled by harmony [not!]

**R – Xi'an's Og**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

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

**I**n 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

**leave a comment**for the author, please follow the link and comment on their blog:

**R – Xi'an's Og**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

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