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()