The below python code implements the Metropolis algorithm and samples from a single variable gaussian distribution. The initial value is sampled uniformly within 5 standard deviations of the mean. Following perturbations are sampled uniformly (+/- 1 standard deviation) and added to current value. A random event is generated in range [0,1], if this value is less than the likelihood ratio of proposed/current, the movement is executed. Otherwise, the current is maintained for another iteration.
Because I'm sampling perturbations from a uniform distribution, inherently, symmetric, I'm just executing the Metropolis algorithm. I'd like to understand MH better, which makes use of (and accounts for) non-symmetric proposal distributions. A few questions:
(1) Why would we want to sample from a non-symmetric proposal distribution and can you provide a concrete example of one (which would take the place of the random.uniform(0,1)
line)?
(2) Can you alter the code detailed below to change M -> MH, and make use of the proposal distribution in the answer to (1) above?
thank you!
def normal(x,mu,sigma):
numerator = np.exp((-(x-mu)**2)/(2*sigma**2))
denominator = sigma * np.sqrt(2*np.pi)
return numerator/denominator
def gaussian_mcmc(hops,mu,sigma):
states = []
burn_in = int(hops*0.2)
current = random.uniform(-5*sigma+mu,5*sigma+mu)
for i in range(hops):
states.append(current)
movement = current + random.uniform(-1,1)
curr_prob = normal(x=current,mu=mu,sigma=sigma)
move_prob = normal(x=movement,mu=mu,sigma=sigma)
acceptance = move_prob/curr_prob
event = random.uniform(0,1)
if acceptance > event:
current = movement
return states[burn_in:]
dist = gaussian_mcmc(100_000,mu=0,sigma=1)
plt.hist(dist,normed=1,bins=20)
plt.plot(lines,normal_curve)