将测量误差添加到 pymc 模型答案

【问题标题】：Adding measurement errors to pymc model将测量误差添加到 pymc 模型
【发布时间】：2014-03-07 08:12:28
【问题描述】：

我在pymc2中有以下模型：

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，我可以从我的模拟数据 sample 中恢复参数 alpha 和 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。

非常感谢任何有关制作这项工作的指导。

【问题讨论】：

标签： pymc

【解决方案1】：

看来改变就够了

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

由

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

采样数据是dist 的合理初始值。无论如何，我不明白其中的逻辑，因为其他初始值（例如零数组）带回了无法再次采样alpha 和beta 的问题。

【讨论】：