2

我正在尝试创建与第 175 页 http://www.people.fas.harvard.edu/~plam/teaching/methods/mcmc/mcmc.pdf相同的 Gibbs 采样器， 它是用 R 编写的，但我'正在尝试在 python 中做到这一点。

我的代码是

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

但是，我的样本很快就会变成无穷大，这很糟糕。谁能看到我在这里出错的地方？我是否以正确的方式从 Gamma 分布中采样？

谢谢

4

1 回答 1

6

问题是 numpy.random.gamma 使用与 R 的 rgamma 函数默认不同的 gamma 分布参数化。numpy.random.gamma 的参数是形状和比例，而 rgamma 函数可以采用形状和速率（它也可以采用比例，但您的代码正在使用速率）。您可以通过反转速率将其转化为比例。这是固定代码：

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()
于 2013-07-23T17:29:21.050 回答