# Proposal function for MCMC - prior distribution NORMAL

i have a prior distribution ~N (0, 0.0168) and wanted to use MCMC to sample from the random generations. I have been having trouble interpreting the concept of a proposal function. To be precise if i use a proposal function with N(0 , 1) i get a sample that's too far off the prior and using something closer to the prior like N( 0 , 0.01) or N(0, 0.02) then the results are closer to the target distribution.

But isn't the whole point of a proposal function to generate samples closer to target distribution so that i could ascertain mean and SD from it to get an idea about the future distribution ?

I am using the HASTINGS algorithm and a uniform distribution [0, 1] to generate the random acceptance criteria (and using code samples from here http://mlwhiz.com/blog/2015/08/19/MCMC_Algorithms_Beta_Distribution/)

    import random
import math
import IPython
import numpy as np
import pylab as pl
import scipy.special as ss
import matplotlib
from scipy.stats import norm
pl.rcParams['figure.figsize'] = (10.0, 10.0)

### target distribution for now would be a straightforward Normal distribution with mean 0 and sigma 0.0158

target_mean = 0
target_sigma = 0.0158

### the mean below might not be used but sigma would since the proposal function is normal of the form N( x(t-1), proposal_sigma)
proposal_mean = 0
proposal_sigma = 0.2

def target_distr_s(x):
return (1/ ((target_sigma)*np.sqrt(np.pi)) )*math.exp( -((x - target_mean)**2)/(2*(target_sigma**2)) )

# This Function returns True if the coin with probability P of heads comes heads when flipped.
def random_coin(p):
unif = random.uniform(0,1)
if unif>=p:
return False
else:
return True

# This Function runs the MCMC chain for Beta Distribution.
def norm_mcmc(N_hops):
states = []
### start the current state with some initialization (RANDOM from N(0,1))
current_state = np.random.normal( proposal_mean, proposal_sigma)
print current_state
for i in range(0,N_hops):
states.append(current_state)
### the next state will always depend on just the prior state. So if it hasn't changed the norm function will return
### yet another random number around the same mean and sigma but the target distribution will decide if the density of
### of this number is high enough for the state to be accepted or not

next_state = np.random.normal( current_state ,proposal_sigma)
#        print ('IN NORM_MCMC '+ next_state + ' '+ current_state+' ret value '+target_distr_s( next_state )+' curr state '+target_distr_s( current_state ) )
if(target_distr_s( current_state )>0):
#            print ('HOLA! '+target_distr_s( next_state )/target_distr_s( current_state ))
ap = min( (target_distr_s( next_state )/target_distr_s( current_state )) ,1) # Calculate the acceptance probability

if ap > 1:
current_state = next_state

return states[-1000:] # Returns the last 100 states of the chain

### call the main functions now

def target_dist(x):
return (1/ ((target_sigma)*np.sqrt(np.pi)) )*math.exp( -((x - target_mean)**2)/(2*(target_sigma**2)) )

def plot_exp():
Ly = []
Lx = []
i_list = np.mgrid[-1:1:100j]
for i in i_list:
Lx.append(i)
Ly.append( target_dist (i))
pl.plot(Lx, Ly, label="Real Distribution: ")
pl.hist(norm_mcmc(10000),normed=True,bins =15, histtype='bar',label="Simulated_MCMC:")
pl.legend()
pl.show()

plot_exp()


As far as I can see from your code, you may only be accepting a proposal when the target density is greater at the proposal value than at the previous value (if ap > 1: current_state = next_state), but you should also accept with probability ap when 0<=ap<=1.

When you use a proposal function that is not optimal (in this case the standard deviation for a N(0,1) distribution is much larger than for the target distribution so that you would get a low acceptance rate of moves), this may show up in MCMC diagnostics (e.g. trace plots, autocorrelation plots, effective sample size, Geweke diagnostic, Gelman-Rubin diagnostic etc.). Did you check any of those? A simple plot of the value of your chain over time (trace plot) could be very helpful. There are some ways of tuning the proposal distribution during an initial stage, but in simple cases like this it might be enough to just generate a longer chain and to only retain every 10th or 100th or 1000th sample.

• Thanks a ton Bjorn..i had forgotten to add the check for accepting with ap if it were in between 0 and 1 ... but before that i went ahead and changed ap>1 to ap>=1 and it seemed to match the target distribution in most plots i took ..so adding the "else accept if its met with prob ap" didn't help much ..any idea why this would be so ? so my mcmc sampler works fine now but its mostly due to changing > 1 to >=1 ..only i am unable to grasp the significance..pointers would be greatly appreciated. – Vikram Murthy Sep 28 '15 at 11:42
• I cannot see how a sampler that only ever moves to higher density values would ever work - would it not just converge towards the mode of the density? Or is there some rounding going on??? Did you try diagnostic plots to see what is going on? – Björn Sep 28 '15 at 12:31

As per your calculation of $ap$(acceptance probability) it would never be greater than 1. And you check for $ap>1$ here.

You are doing ap = min( (target_distr_s( next_state )/target_distr_s( current_state )) ,1) so ap would be <=1

You should flip a coin with probability $ap$ and move to the next state only if coin lands heads.

• Thanks for the helpful post and the feedback Rahul , but according to the code if the probability distribution for the next random variable is greater than the probability distribution for the current one, surely the ratio would be greater than 1 , hence making the check work no ? i say this because though my code was wrong in that i forgot to add the checks w.r.t. "accepting with prob ap" and using just >1 instead of >=1 as specified in the metropolis hastings algorithm, my plots were showing some ok results .. – Vikram Murthy Sep 28 '15 at 11:46
• See how we calculate ap here: ap = min( (target_distr_s( next_state )/target_distr_s( current_state )) ,1) . ap is a probability and so i made it a point of keeping it between 0 and 1 and hence i use the min function which makes sure that 0<=ap<=1. So if the probability distribution for the next random variable is greater than the probability distribution for the current one, surely the ratio would be greater than 1, and the min function would make ap=1. – Rahul Agarwal Sep 28 '15 at 13:29
• just change [if ap > 1:] line to [ if random_coin(ap):] – Rahul Agarwal Sep 28 '15 at 13:35