1

同じインターセプトを共有する複数の行を合わせようとしています。

import numpy as np
import pymc

# Observations
a_actual = np.array([[2., 5., 7.]]).T
b_actual = 3.
t = np.arange(100)
obs = np.random.normal(a_actual * t + b_actual)


# PyMC Model
def model_linear():
    b = pymc.Uniform('b', value=1., lower=0, upper=200)

    a = []
    s = []
    r = []
    for i in range(len(a_actual)):
        s.append(pymc.Uniform('sigma_{}'.format(i), value=1., lower=0, upper=100))
        a.append(pymc.Uniform('a_{}'.format(i), value=1., lower=0, upper=200))
        r.append(pymc.Normal('r_{}'.format(i), mu=a[i] * t + b, tau=1/s[i]**2, value=obs[i], observed=True))

    return [pymc.Container(a), b, pymc.Container(s), pymc.Container(r)]

model = pymc.Model(model_linear())
map = pymc.MAP(model)
map.fit()
map.revert_to_max()

計算された MAP 推定値は、実際の値とはかけ離れています。これらの値は、 と の下限と上限、 の実際の値 (たとえば、適切な見積もりが得られる)、または回帰を実行する行数にも非常sigmasにa敏感aですa = [.2, .5, .7]。

これは私の線形回帰を実行する正しい方法ですか?

ps : シグマに指数事前分布を使用しようとしましたが、結果は良くありませんでした。

4

1 に答える 1

1

MAP を使用するのは最善の策ではないかもしれません。適切なサンプリングを行うことができる場合は、サンプル コードの最後の 3 行を次のように置き換えることを検討してください。

MCMClinear = pymc.MCMC( model)
MCMClinear.sample(10000,burn=5000,thin=5)
linear_output=MCMClinear.stats()

これを印刷するlinear_outputと、パラメーターの非常に正確な推論が得られます。

于 2014-04-10T12:37:56.517 に答える