Python吉布斯採樣器不工作

Question

我試圖創建相同的吉布斯採樣器，在這裏找到175頁 http://www.people.fas.harvard.edu/~plam/teaching/methods/mcmc/mcmc.pdf 這是寫在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分佈中抽樣？

感謝

Answer 1

的問題是，numpy.random.gamma使用除R的rgamma功能默認完成伽瑪分佈的不同參數。 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()