PyMC多重線性迴歸

我想要適合幾條線共享相同的截距。PyMC多重線性迴歸

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：我試圖用指數先驗分佈來計算sigma，但結果並不好。

來源

2014-04-07 synapski

你不能只運行3個單獨的線性迴歸？ – nstjhp

我不能每個觀察都需要共享相同的截距（'b'）。我的真實模型比這更復雜;這是我的問題的簡化版本。 – synapski

我認爲使用MAP可能不是您最好的選擇。如果你能夠做一個適當的採樣，然後考慮

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

打印linear_output此給出的參數非常準確的推論將最後3行的示例代碼。

來源

2014-04-10 12:37:56 nstjhp

是的，MCMC採樣工作正常，效果良好。我只是想知道爲什麼MAP估計值對變量界限或對簡單模型觀測值的實際值非常敏感。使用[lmfit]（http://cars9.uchicago.edu/software/python/lmfit/）庫（它也基於scipy的優化包）對同樣的問題給出了很好的結果。 – synapski

PyMC多重線性迴歸

回答

相關問題