添加測量誤差pymc模型

我在pymc2以下模型：添加測量誤差pymc模型

import pymc 
from scipy.stats import gamma 

alpha = pymc.Uniform('alpha', 0.01, 2.0) 
scale = pymc.Uniform('scale', 1.0, 4.0) 

@pymc.deterministic(plot=False) 
def beta(scale=scale): 
    return 1.0/scale 

@pymc.potential 
def p_factor(alpha=alpha, scale=scale, lmin=lmin, n=len(sample)): 
    dist = gamma(alpha, loc=0., scale=scale) 
    fp = 1.0 - dist.cdf(lmin) 
    return -(n+1)*np.log(fp) 

obs = pymc.Gamma("obs", alpha=alpha, beta=beta, value=sample, observed=True)

該模型的物理背景是luminosity function of galaxies（LF），即，具有亮度L.對於一些星系的概率LF只是一種伽馬函數。由於星系調查通常會錯過很大一部分目標，特別是那些低亮度的目標，因此可能會造成數據截斷。在此模型中，我想念下面的所有內容lmin

此方法的詳細信息可在this paper by Kelly et al中找到。

這個模型的工作：我在模型上運行MAP和MCMC，我可以恢復的參數從我的模擬數據samplealpha和scale，增加不確定性lmin增長。

現在我想插入高斯測量誤差。爲了簡單起見，所有數據具有相同的精度。我也沒有修改包括錯誤的可能性。

alpha = pymc.Uniform('alpha', 0.01, 2.0) 
scale = pymc.Uniform('scale',1.0, 4.0) 
sig = 0.1 
tau = math.pow(sig, -2.0) 

@pymc.deterministic(plot=False) 
def beta(scale=scale): 
    return 1.0/scale 

@pymc.potential 
def p_factor(alpha=alpha, scale=scale, lmin=lmin, n=len(sample)): 
    dist = gamma(alpha, loc=0., scale=scale) 
    fp = 1.0 - dist.cdf(lmin) 
    return -(n+1) * np.log(fp) 

dist = pymc.Gamma("dist", alpha=alpha, beta=beta) 
obs = pymc.Normal("obs", mu=dist, tau=tau, value=sample, observed=True)

但可以肯定，我在這裏做得不對，因爲這種模式是行不通的。當我在這個模型上運行pymc.MAP我收回alpha和初始值scale

vals = {'alpha': alpha, 'scale': scale, 'beta': beta, 
    'p_factor': p_factor, 'obs': obs, 'dist': dist} 
M2 = pymc.MAP(vals) 
M2.fit() 
print M2.alpha.value, M2.scale.value 
>>> (array(0.010000000006018368), array(1.000000000833973))

當我運行pymc.MCMC，alpha和beta沒有跟蹤的。

M = pymc.MCMC(vals) 
M.sample(10000, burn=5000) 
... 
M.stats()['alpha'] 
>>> {'95% HPD interval': array([ 0.01000001, 0.01000502]), 
'mc error': 2.1442678276712383e-07, 
'mean': 0.010001588137798096, 
'n': 5000, 
'quantiles': {2.5: 0.0100000088679046, 
25: 0.010000382359859467, 
50: 0.010001100377476166, 
75: 0.010001668672799679, 
97.5: 0.0100050194240779}, 
'standard deviation': 2.189828287191421e-06}

再次初始值。實際上，如果我將alpha更改爲0.02，那麼alpha的恢復值爲0.02。

這是a notebook with the working model plus simulated data。

這是a notebook with the error model plus simulated data。

任何關於使這項工作的指導將非常感激。

來源

2014-02-25 Sergio

似乎是足以改變

dist = pymc.Gamma("dist", alpha=alpha, beta=beta)

通過

dist = pymc.Gamma("dist", alpha=alpha, beta=beta, value=sample)

的採樣數據是用於dist合理的初始值。無論如何，我沒有得到邏輯，因爲其他初始值（例如一個零數組）帶回了再次不採樣alpha和beta的問題。

來源

2014-03-07 00:06:51 Sergio

添加測量誤差pymc模型

回答

相關問題