0

私は長い間 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 の線形回帰 (被験者ごとに 1 つ) を設定し、切片パラメーターのベクトルを最終的な回帰のデータとして使用することでした。

#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)

このモデルは、サンプリング ステップで失敗します。dtype=object の代わりに、offsetArr を dtype=float64 としてキャストしようとすると、エラーが発生するようです。これを行う正しい方法は何ですか？私のoffsetArrをfloat64にキャストするために別の決定論的ノードが必要ですか? 特別な種類の pymc.Container が必要ですか? ご協力いただきありがとうございます！

4

1 に答える 1

0

PyMC オブジェクトを格納するために、numpy 配列の代わりに単純なリストを使用してみましたか?

于 2013-06-11T16:02:44.963 に答える