许多回归的PyMC回归?(PyMC regression of many regressions?)

2019-10-29 12:30发布

我没有使用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的? 谢谢你的帮助!

Answer 1:

你是否尝试过使用一个简单的列表,而不是numpy的阵列来存储PyMC对象?



文章来源: PyMC regression of many regressions?