PyMC 多重线性回归

Question

我正在尝试拟合几条共享相同截距的直线。

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

计算出来的最优估计值和实际值差得很远。而且这些值对sigmas和a的上下限非常敏感，a的实际值也会影响结果（比如说a = [.2, .5, .7]会给我比较好的估计），或者说回归的直线数量也会有影响。

这样做线性回归的方法对吗？

补充：我尝试使用指数先验分布来处理sigmas，但结果并没有更好。

Answer 1

我觉得用MAP可能不是最好的选择。如果你能进行合适的抽样，那就考虑把你示例代码的最后三行换成下面的内容：

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

打印这个linear_output可以得到非常准确的参数推断。