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嗎?謝謝你的幫助!