PYMC2 零概率错误

PYMC2 ZeroProbability error

我有一个具有 6 个参数且具有统一先验的模型:

parameter1 = pm.Uniform('parameter1',0.01,1)
parameter2 = pm.Uniform('parameter2',0,2)
parameter3 = pm.DiscreteUniform('parameter3',1,50)
parameter4 = pm.Uniform('parameter4',0,1.75)
parameter5 = pm.Uniform('parameter5', 0.005, 0.25)
parameter6 = pm.Uniform('parameter6', 0.005, 0.15)

我有一个自定义似然函数 returns 对数似然值:

@pm.potential
def log_l(experiment=experiment,parameter1=parameter1,parameter2=parameter2,parameter3=parameter3,parameter4=parameter4,parameter5=parameter5,parameter6=parameter6):

    if parameter5<parameter4:
        return -np.inf

    parameters=[parameter1, parameter2, parameter3]

    log_l=calculate_probability(parameters, t_m, tol, n_integration, parameter4, parameter5, parameter6, experiment.decon_all[freq,:,:])

    return log_l

其中 calculate_probability 是我的函数,它 returns 在给定参数值和观测数据的情况下记录该模型的可能性。 由于某些原因,当 MCMC 采样时:

model = pm.MCMC([parameter1,parameter2,parameter3,parameter4,parameter5,parameter6,log_l])
model.sample(100)

并且程序满足 if 条件 (parameter5<parameter4) 我得到这个错误:

pymc.Node.ZeroProbability: Potential log_l forbids its parents' current values

我想知道是否有人知道我可能做错了什么?

根据 the docs:

Potentials have one important attribute, logp, the log of their current probability or probability density value given the values of their parents. The only other additional attribute of interest is parents, a dictionary containing the potential’s parents.

所以你的 log_l 的 parent 是 log_l 的参数:

In [11]: list(log_l.parents.keys())
Out[11]: ['experiment', 'parameter2', 'parameter3', 'parameter1', 'parameter4', 
          'parameter6', 'parameter5']

根据 this answer(我强调):

When a random variable is defined as a function of another random variable, PyMC checks that no value of the parent distribution leads to an impossible value for the child distribution.

使log_lreturn-np.inf意味着parent变量的某些值是不可能的。因此,PyMC 引发了一个 ZeroProbability 异常。


所以不要使用

来约束模型
if parameter5<parameter4:
    return -np.inf

你可以用

定义parameter5
parameter4 = pm.Uniform('parameter4', 0, 1.75)
parameter5 = pm.Uniform('parameter5', parameter4, 0.25)

确保parameter5 > parameter4.