<p><a target="_blank" rel="nofollow" href="https://i0.wp.com/xianblog.files.wordpress.com/2012/08/img_1974.jpg?ssl=1"><img sizes="(max-width: 450px) 100vw, 450px" data-image-description="" data-orig-file="https://i0.wp.com/xianblog.files.wordpress.com/2012/08/img_1974.jpg?resize=450%2C338&amp;ssl=1" data-orig-size="2816,2112" data-image-title="(past) sunset over Kata Tjuta (Mounts Olga); Northern Territory, Australia, Aug. 6, 2012" class="wp-image-17249" data-comments-opened="1" data-recalc-dims="1" data-attachment-id="17249" title="(past) sunset over Kata Tjuta (Mounts Olga); Northern Territory, Australia, Aug. 6, 2012" data-permalink="https://xianblog.wordpress.com/2012/08/12/kata-tjuta-sunset/img_1974/" data-medium-file="https://i0.wp.com/xianblog.files.wordpress.com/2012/08/img_1974.jpg?resize=450%2C338&amp;ssl=1?w=300" data-large-file="https://i0.wp.com/xianblog.files.wordpress.com/2012/08/img_1974.jpg?resize=450%2C338&amp;ssl=1?w=450" data-image-meta="{&quot;aperture&quot;:&quot;3.5&quot;,&quot;credit&quot;:&quot;&quot;,&quot;camera&quot;:&quot;Canon IXUS 115 HS&quot;,&quot;caption&quot;:&quot;&quot;,&quot;created_timestamp&quot;:&quot;1344251437&quot;,&quot;copyright&quot;:&quot;&quot;,&quot;focal_length&quot;:&quot;8.022&quot;,&quot;iso&quot;:&quot;640&quot;,&quot;shutter_speed&quot;:&quot;0.0125&quot;,&quot;title&quot;:&quot;&quot;}" src="https://i0.wp.com/xianblog.files.wordpress.com/2012/08/img_1974.jpg?resize=450%2C338&ssl=1"></a> <strong>I</strong>n answering <a target="_blank" rel="nofollow" href="http://stats.stackexchange.com/a/248697/7224">yet another question</a> on X validated about the numerical approximation of the marginal likelihood, I suggested using an <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/2016/07/15/the-curious-incident-of-the-inverse-of-the-mean/">harmonic mean estimate</a> 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.</p>
<pre><code>f &lt;- function(x){
    if ((0&lt;x)&amp;(x&lt;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) &lt; aprob){
        x = can
        fx = fcan}
    chain=c(chain,fx)}
I = pi*length(chain)/sum(1/chain) #integral harmonic approximation</code></pre>
<p>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!</p>
<p>Filed under: <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/books/">Books</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/kids/">Kids</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/mountains/">Mountains</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/pictures/">pictures</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/statistics/r-statistics/">R</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/running/">Running</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/statistics/">Statistics</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/category/travel/">Travel</a> Tagged: <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/gaussian-random-walk/">Gaussian random walk</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/harmonic-mean-estimator/">harmonic mean estimator</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/metropolis-hastings-algorithm/">Metropolis-Hastings algorithm</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/monte-carlo-statistical-methods/">Monte Carlo Statistical Methods</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/numerical-integration/">numerical integration</a>, <a target="_blank" rel="nofollow" href="https://xianblog.wordpress.com/tag/simulation/">simulation</a> <a target="_blank" rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/xianblog.wordpress.com/32177/"><img alt="" src="http://feeds.wordpress.com/1.0/comments/xianblog.wordpress.com/32177/" border="0"></a> <img alt="" src="https://i1.wp.com/pixel.wp.com/b.gif?resize=1%2C1&ssl=1" border="0" data-recalc-dims="1"></p>
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/Bi0GsqvXQek">