How do I sample from a Gamma(26,6) distribution using Metropolis Hastings

user_11628437 0
Viewed 101

I am studying MCMC and decided to test my understanding of Metropolis Hastings algorithm on a conjugate Bayesian problem. Given a poisson likelihood and gamma prior I am trying to plot the posterior distribution. The data is given as follows, y = (4, 4, 5, 8, 3) and the gamma prior parameters are a = b = 1.

Since, this is a conjugate model, I came up with the true posterior, Ga(a + n*y_bar, b + n) = Ga(25, 6). Here is the code for it -

##### The real distribution ####
import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
data = np.array([4, 4, 5, 8, 3])
n = len(data)
#Gamma(shape = 25, rate = 6)
lik = st.gamma.pdf(a = 25, scale = 1/6, x = np.linspace(0, 10, num=100))
plt.plot(lik)
plt.show()

Here is my Metropolis Hastings code. However, I am unable to get the true posterior using it. I borrowed a lot of ideas from here - https://people.duke.edu/~ccc14/sta-663/MCMC.html

def target(lambda_val, log_val = True):
    if log_val == False:
        return ((lambda_val)**(n*np.mean(data)+a+1))*np.exp(-lambda_val*(b+n))
    else:
        return np.exp((n*np.mean(data)+a+1)*np.log(lambda_val) - lambda_val*(n+b))

### Using MCMC to sample from it #####
niters = 1000
theta = 40

a = 1
b = 1
samples = np.zeros(niters+1)
samples[0] = 5
naccept = 0
for i in range(niters):
    theta_p = theta + st.norm(0,3).rvs()
    rho= min(1, target(theta_p)/target(theta))
    u = np.random.uniform()
    if u<rho:
        naccept += 1
        theta = theta_p
    samples[i+1] = theta
nmcmc = len(samples)//2
print("Efficiency = ", naccept/niters)

plt.hist(samples[nmcmc:], 40)
plt.show()

Please let me know if I missed any details.

mcmc bayesian python
1 Answers

The likelihood is not properly defined. Specifically, one needs to explicitly exclude proposals outside the support (non-positives) by returning 0 (or -Inf for logspace).

Also the step size (sd=3) is rather large considering the actual location (~4). This results in frequently proposing negative numbers.

Lastly, why is the initialization theta=40? Perhaps this resulted from misinterpreting the plot, which doesn't use a proper x-axis, but only an index:

enter image description here

Plotting with an valid x-axis:

xs = np.linspace(0, 10, 100)
lik = st.gamma.pdf(a=25, scale=1/6, x=xs)
plt.plot(xs, lik)
plt.show()

enter image description here

one sees that a more appropriate initialization might be 4.

user_5709180
Related