我使用 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 转换为 dtype=float64 而不是它的 dtype=object 时。这样做的正确方法是什么？我是否需要另一个确定性节点来帮助将我的 offsetArr 转换为 float64？我需要一种特殊的 pymc.Container 吗？谢谢你的帮助！