2013-05-07 88 views
0

我還沒有被使用PyMC時間不長,但我很高興在我的速度有多快能夠功成線性迴歸(此代碼應該沒有IPython的修改即可運行):許多回歸的PyMC迴歸?

import pandas as pd 
from numpy import * 
import pymc 

data=pd.DataFrame(rand(40)) 
predictors=pd.DataFrame(rand(40,5)) 
sigma = pymc.Uniform('sigma', 0.0, 200.0, value=20) 
params= array([pymc.Normal('%s_coef' % (c), mu=0, tau=1e-3,value=0) for c in predictors.columns]) 

@pymc.deterministic(plot=False) 
def linear_regression_model(x=predictors,beta=params): 
    return dot(x,beta) 

ynode = pymc.Normal('ynode', mu=linear_regression_model, tau=1.0/sigma**2,value=data,observed=True) 


m = pymc.Model(concatenate((params,array([sigma,ynode])))) 

%time pymc.MCMC(m).sample(10000, 5000, 5, progress_bar=True) 

在這個模型中,每個科目有40個科目(觀察)和5個協變量。由於隨機數據,模型不會收斂,但它沒有錯誤的樣本(我的真實數據確實收斂到了精確的結果)。

我遇到問題的模型是此的擴展。每個主題實際上有3個(或N個)觀測值,所以我需要對這些觀測值擬合一條線,然後使用該線的截距作爲最終迴歸節點的「數據」或輸入。我認爲這是一個經典的等級模型,但如果我以錯誤的方式思考它,則糾正我。

我的解決方案是建立一系列40個線性迴歸(每個主題一個迴歸),然後使用截距參數向量作爲最終迴歸的數據。

#nodes for fitting 3 values for each of 40 subjects with a line 
#40 subjects, 3 data points each 
data=pd.DataFrame(rand(40,3)) 
datax=arange(3) 

""" 
to fit a line to each subject's data we need: 
    (1) a slope and offset parameter 
    (2) a stochastic node for the data 
    (3) a sigma parameter for the stochastic node 

Initialize them all as object arrays 
""" 
sigmaArr=empty((len(data.index)),dtype=object) 
ynodeArr=empty((len(data.index)),dtype=object) 
slopeArr=empty((len(data.index)),dtype=object) 
offsetArr=empty((len(data.index)),dtype=object) 

#Fill in the empty arrays 
for i,ID in enumerate(data.index): 
    sigmaArr[i]=pymc.Uniform('sigma_%s' % (ID) , 0.0, 200.0, value=20) 
    slopeArr[i]=pymc.Normal('%s_slope' % (ID), mu=0, tau=1e-3,value=0) 
    offsetArr[i]=pymc.Normal('%s_avg' % (ID), mu=0, tau=1e-3,value=data.ix[ID].mean()) 

    #each model fits a line to the three data points 
    @pymc.deterministic(name='time_model_%s' % ID,plot=False) 
    def line_model(xx=datax,slope=slopeArr[i],avg=offsetArr[i]): 
     return slope*xx + avg 

    ynodeArr[i]=pymc.Normal('ynode_%s' % (ID), mu = line_model, tau = 1/sigmaArr[i]**2,value=data.ix[ID],observed=True) 


#nodes for final regression (there are 5 covariates in this regression) 
predictors=pd.DataFrame(rand(40,5)) 
sigma = pymc.Uniform('sigma', 0.0, 200.0, value=20) 
params= array([pymc.Normal('%s_coef' % (c), mu=0, tau=1e-3,value=0) for c in predictors.columns]) 


@pymc.deterministic(plot=False) 
def linear_regression_model(x=predictors,beta=params): 
    return dot(x,beta) 

ynode = pymc.Normal('ynode', mu=linear_regression_model, tau=1.0/sigma**2,value=offsetArr) 

nodes=concatenate((sigmaArr,ynodeArr,slopeArr,offsetArr,params,array([sigma, ynode]))) 

m = pymc.Model(nodes) 

%time pymc.MCMC(m).sample(10000, 5000, 5, progress_bar=True) 

該模型在採樣步驟失敗。該錯誤似乎是在試圖將offsetArr轉換爲dtype = float64而不是其dtype = object時。什麼是正確的方法來做到這一點?我是否需要另一個確定性節點來幫助將offsetArr轉換爲float64?我需要一種特殊的pymc.Container嗎?謝謝你的幫助!

回答

0

您是否嘗試過使用簡單列表而不是numpy數組來存儲PyMC對象?