データから単純な線形関数とガンマ分布ノイズ項のパラメーターを推定したいと思います。(注: これはhttps://stats.stackexchange.com/questions/88676/regression-with-unidirectional-noiseのフォローアップの質問ですが、単純化され、より実装固有です)。次のように生成された観測データがあるとします。
import numpy as np
np.random.seed(0)
size = 200
true_intercept = 1
true_slope = 2
# Generate observed data
x_ = np.linspace(0, 1, size)
true_regression_line = true_intercept + true_slope * x_ # y = a + b*x
noise_ = np.random.gamma(shape=1.0, scale=1.0, size=size)
y_ = true_regression_line + noise_
次のようになります。
次のように、pymc を使用してこれらのパラメーターを推定してみました。
from pymc import Normal, Gamma, Uniform, Model, MAP
# Define priors
intercept = Normal('intercept', 0, tau=0.1)
slope = Normal('slope', 0, tau=0.1)
alpha = Uniform('alpha', 0, 2)
beta = Uniform('beta', 0, 2)
noise = Gamma('noise', alpha=alpha, beta=beta, size=size)
# Give likelihood > 0 to models where the regression line becomes larger than
# any of the datapoint
y = Normal('y', mu=intercept + slope * x_ + noise, tau=100,
observed=True, value=y_)
# Perform MAP fit of model
model = Model([alpha, beta, intercept, slope, noise])
map_ = MAP(model)
map_.fit()
ただし、これにより、真の値からかけ離れた見積もりが得られます。
- 切片: 真: 1.000、推定: 3.281
- 勾配: 真: 2.000、推定: -3.400
私は何か間違ったことをしていますか？