Python Gibbs Sampler not working

Question

I'm trying to create the same Gibbs sampler that is found here on page 175 http://www.people.fas.harvard.edu/~plam/teaching/methods/mcmc/mcmc.pdf which was written in R but I'm trying to do it in python.

my code is

from numpy import *
import matplotlib.pylab as pl

def gibbs_sampler(alpha,delta,gamma,y,t):
    #initialize beta
    beta=1

    num_iter=100

    beta_draws=[]
    lambda_draws=[]

    for i in range(num_iter):
        #sample lambda given other lambdas and beta
        lambdas=lambda_update(alpha,beta,y,t)

        #record sample
        lambda_draws.append(lambdas)

        #sample beta given lambda samples
        beta=beta_update(alpha,gamma,delta,lambdas,y)

        #record sample
        beta_draws.append(beta)


    pl.plot(array(beta_draws))
    pl.show()

def lambda_update(alpha,beta,y,t):
    new_alpha=[(x+alpha) for x in y]
    new_beta=[(a+beta) for a in t]

    #sample from this distribution 10 times
    samples=random.gamma(new_alpha,new_beta)
    return samples


def beta_update(alpha,gamma,delta,lambdas,y):
    #get sample
    sample=random.gamma(len(y)*alpha+gamma,delta+sum(lambdas))
    return sample



def main():
    y=[5,1,5,14,3,19,1,1,4,22]
    t=[94,16,63,126,5,31,1,1,2,10]


    alpha=1.8
    gamma=0.01
    delta=1

    gibbs_sampler(alpha,delta,gamma,y,t)

if __name__ == '__main__':
    main()

however, my samples quickly go to infinity which is bad. Can anyone see where I am going wrong here? Am I sampling from the Gamma distribution in the correct way?

Thanks

Answer 1

The problem is that numpy.random.gamma uses a different parameterization of the gamma distribution than R's rgamma function does by default. The parameters of numpy.random.gamma are shape and scale, whereas the rgamma function can take shape and rate (it can also take scale, but your code is using the rate). You can turn rate into scale by inverting it. Here is the fixed code:

from numpy import *
import matplotlib.pylab as pl

def gibbs_sampler(alpha,delta,gamma,y,t):
    #initialize beta
    beta=1

    num_iter=100

    beta_draws=[]
    lambda_draws=[]

    for i in range(num_iter):
        #sample lambda given other lambdas and beta
        lambdas=lambda_update(alpha,beta,y,t)

        #record sample
        lambda_draws.append(lambdas)

        #sample beta given lambda samples
        beta=beta_update(alpha,gamma,delta,lambdas,y)

        #record sample
        beta_draws.append(beta)

    pl.plot(beta_draws)
    pl.show()

def lambda_update(alpha,beta,y,t):

    new_alpha=[(x+alpha) for x in y]
    new_beta=[1.0/(a+beta) for a in t]#Changed here

    #sample from this distribution 10 times
    samples=random.gamma(new_alpha,new_beta)
    return samples


def beta_update(alpha,gamma,delta,lambdas,y):
    #get sample
    sample=random.gamma(len(y)*alpha+gamma, 
                        1.0 / (delta+sum(lambdas)))#Changed here
    return sample


def main():

    y=[5,1,5,14,3,19,1,1,4,22]
    t=[94,16,63,126,5,31,1,1,2,10]

    alpha=1.8
    gamma=0.01
    delta=1

    gibbs_sampler(alpha,delta,gamma,y,t)

if __name__ == '__main__':
    main()