我没有使用PyMC时间不长,但我很高兴在我的速度有多快能得到一元线性回归掉在地上(此代码应在不修改的IPython运行):
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为D类= float64时,而不是它的D型=对象。 什么是这样做的正确方法? 我需要另一个确定性节点,以帮助投我offsetArr到float64? 我需要一种特殊的pymc.Container的? 谢谢你的帮助!