3
$\begingroup$

I've been asked to approximate the integral of a one dimensional unnormalised posterior with a flat prior, using a Metropolis Hastings Markov Chain Monte Carlo, I realise that this isn't a practical endeavour, i think its more of an instructional exercise.

At the moment, I am creating a Markov chain of sampling points, $\{x_i\}$, for the distribution using the MH algorithm, and then approximating the integral of the posterior as:

\begin{equation} I \approx (\text{Prior interval}) * \frac{1}{N}\sum_i f(x_i) \end{equation}

where $f(x_i)$ is the likelihood evaluated at the point $x_i$ in the chain, and N is the length of the Markov Chain.

I have tested this algorithm in Python for numerous simple posteriors, and it converges on completely the wrong answer in every case. I've been over the code many times, and I'm almost certain there's nothing wrong with it. I also read that you need to discard the first half of the chain, and discard everything but every, say, 10th point. I've tried this and it makes no difference.

The code is shown below. It tries to integrate $\sin x$ between $0$ and $\pi$. For 20000 iterations it gives $2.45448671306 \pm 0.00502441977146$, which is obviously very wrong (answer should be 2) .The graph below shows the convergence of the algorithm convergence graph

def f(x):
    if 0. <= x <= np.pi:
        return np.sin(x)
    else:
        return 0


n = 20000 #number of iterations
sigma = 0.5

x = uniform(0, np.pi) #initial x value

chain = [] #stores Markov chain 
chain.append(f(x)) #initialises chain



accept = 0 #number of accepted jumps

#generates an array of random x values from norm distribution
rands = normal(0, sigma, n) 

#Metropolis - Hastings algorithm
for i in range(1,n):
    can = x + rands[i]  #candidate for jump
    aprob = min([1., f(can)/f(x)]) #acceptance probability
    u = uniform(0,1) #rand number between 0 and 1
    if u < aprob:
        x = can
        accept += 1
    chain.append(f(x))

I = (np.pi)*sum(chain)*(1/len(chain)) #integral approximation

#uncertainty on approx.
uncertainty = (np.pi) / np.sqrt(len(chain)) * np.std(chain) 

#integral approx at each iteration
conv = np.array([(np.pi)*sum(chain[:x+1])*(1/(x+1)) for x in range(len(chain))])

print('MCMC approx: ', I, '\t','uncertainty: ', uncertainty)
$\endgroup$
  • $\begingroup$ As discussed in this answer, using the posterior sample to approximate the marginal likelihood, ie the integrated unnormalised posterior, is not correct. You either need to use a prior sample or the posterior sample with a more convoluted estimate. I list a collection of solutions in this answer. $\endgroup$ – Xi'an Nov 29 '16 at 19:52
  • 1
    $\begingroup$ Thank you for your reply, that answer works like magic! I'll have to go away and read up on why that works. However, the way the question was set, and from my lecture notes, suggests that I really should be able to tackle this problem just using the MH algorithm and the mean of the Markov Chain (we haven't covered anything more sophisticated than this, and aren't expected to). So are you sure that a simpler alteration to the algorithm can't fix it? And could you point me to some literature which explains why the mean of the markov chain wont give the integral? $\endgroup$ – oweydd Nov 29 '16 at 22:14
  • 1
    $\begingroup$ Ok, so I've asked my lecturer and we weren't actually required to find the integral, just to sample from the posterior. Nonetheless, i think I've learnt something useful, so thank you. $\endgroup$ – oweydd Nov 30 '16 at 10:11
  • $\begingroup$ There was a major mistake in your code, namely that you did not account for the rejections in the Metropolis-Hasting part. When a proposed value can is rejected the current value x needs to be replicated, as illustrated in my code. $\endgroup$ – Xi'an Nov 30 '16 at 10:19
5
$\begingroup$

"I've been over the code many times, and I'm almost certain there's nothing wrong with it."

There is an issue with the code, actually, in that you fail to replicate the current value of the Markov chain x when rejecting the candidate value can, while you should: thus the conditional update

   if u < aprob:
        x = can
        accept += 1
    chain.append(f(x))

should be corrected as done in the R code below:

    if (u < aprob){
        x = can
        accept = accept+1}
    chain=c(chain,f(x))}

Furthermore, there is a conceptual (and common) error in averaging the $f(x_i)$'s when the $x_i$'s are simulated from the posterior distribution, $g(x)\propto f(x)\mathbb{I}_{(0,\pi)}(x)$. Indeed, the expectation of a term in this average is$$\int_0^\pi f(x)g(x)\,\text{d}x=\int_0^\pi f(x)\times\frac{f(x)\mathbb{I}_{(0,\pi)}(x)}{\int _0^\pi f(x)\,\text{d}x}\,\text{d}x=\frac{\int_0^\pi f^2(x) \,\text{d}x}{\int _0^\pi f(x)\,\text{d}x}$$

Using the posterior sample as you do in the Monte Carlo approximation hence amounts to using the likelihood twice, i.e. integrating $f^2$ rather than $f$.


The Harmonic Mean Estimator

One immediate way to deduce an approximation of the integral$$\mathcal{I}=\dfrac{1}{\pi}\int_0^\pi f(x)\text{d}x=\dfrac{\cos(\pi)-\cos(0)}{\pi}=\frac{2}{\pi}$$ which is the marginal likelihood, is to use the sample $(x_i)$ from the posterior inside an harmonic mean estimate: $$\hat{\mathcal{I}}=1\big/\frac{1}{n}\sum_{i=1}^n \frac{1}{f(x_i)}$$This is a consequence of the identity$$\mathbb{E}[f(X)^{-1}]=\int_0^\pi f(x)^{-1}\frac{f(x)\mathbb{I}_{(0,\pi}(x)}{\mathcal{I}\pi}\text{d}x=\int_0^\pi\frac{\mathbb{I}_{(0,\pi}(x)}{\mathcal{I}\pi}\text{d}x=\frac{1}{\mathcal{I}}$$ The estimator returns an unbiased estimate of the inverse integral, but its major drawback is that it often has an infinite variance, hence is unreliable no matter what the number of simulations. Radford Neal has called this approach the worst Monte Carlo method ever and I concur. Only use it when you want to check that the value of the integral is what another approach tells you it should be. As in this exercise.

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

Here is the corresponding result for n=1e5:

> I
[1] 2.000591

which coincides with the true value. Hence the harmonic mean did not hit too small a value of $\sin(x)$ to induce large variations in the average. But, nonetheless, $1/\sin(x)$ has an infinite integral on $(0,\pi)$ and the resulting harmonic mean estimator has an infinite variance.

An explanation for the surprisingly good performances of this "worst Monte Carlo method ever" is that sigma = 0.5 does not often return values of $\theta$ that are close to $0$ or to $\pi$. If instead one uses sigma = 0.01 in the above code, the result is quite different:

> I
[1] 3.039536

and highly variable

> I
[1] 0.5901059

Hence it cannot be trusted and should not be used.

$\endgroup$
0
$\begingroup$

The Metro-Hastings sampling procedure is over-sampling higher values due to the accept/reject criteria.

The average of the sample values returned by
sum(chain) * (1/len(chain))
is too high.

In typical method for Monte Carlo integration, random samples are drawn uniformly across the [0, pi] interval.

You can still use a Markov Chain to create "jumps". For example, changing

for i in range(1,n):
can = x + rands[i]  #candidate for jump
aprob = min([1., f(can)/f(x)]) #acceptance probability
u = uniform(0,1) #rand number between 0 and 1
if u < aprob:
    x = can
    accept += 1
chain.append(f(x))

to

for i in range(1, n):
x += rands[i]  #candidate for jump
if f(x) > 0:
    chain.append(f(x))

will get a much better approximation

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.