0

我在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)

该模型的物理背景是星系的光度函数（LF），即星系具有光度L的概率。对于某些类型的星系，LF只是一个伽马函数。可能会导致数据截断，因为星系调查通常会错过很大一部分目标，尤其是那些低光度的目标。在这个模型中，我想念下面的一切lmin

这种方法的细节可以在Kelly 等人的这篇论文中找到。

这个模型有效：我在模型上运行MAP，我可以从我的模拟数据MCMC中恢复参数alpha，随着不确定性的增加。scalesamplelmin

现在我想插入高斯测量误差。为简单起见，所有数据都具有相同的精度。我也没有修改包含错误的可能性。

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的初始值alphascale

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 开始，则恢复值为alpha0.02。

这是一个带有工作模型和模拟数据的笔记本。

这是一个带有误差模型和模拟数据的笔记本。

任何有关使这项工作的指导将不胜感激。

4

1 回答 1

1

似乎足以改变

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

经过

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

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

于 2014-03-07T00:06:51.803 回答