pymc:根据可观察对象的函数推断参数

pymc: Inferring parameters based on functions of observables

我观察到几条光发射线,我有一个模型可以根据我想要的两个参数 qz 预测这些线的几个(通量)比率推断。

我创建了 @pymc.deterministic 个对象,这些对象采用 qz 的值(每个对象在一些物理上感兴趣的区域都没有信息先验),并将它们变成 "predicted"比例。大约有 7 个比率,它们的形式为:

@pymc.deterministic(observed=True, value=NII_SII)
def NII_SII_th(q=q, z=z):
    return NII_SII_g(np.array([q, z]))

我还可以定义从观察中得出的比率,例如

@pymc.deterministic
def NII_SII(NII_6584=NII_6584, SII_6717=SII_6717,
            rcf_NII_6584=rcf_NII_6584, rcf_SII_6717=rcf_SII_6717):
    return np.log10(
        (rcf_NII_6584*NII_6584) / \
        (rcf_SII_6717*SII_6717))

其中,例如,NII_6584 是其中一条线的观测通量,rcf_NII_6584 是同一条线的通量校正。这些校正本身由线波长(以无限精度已知)和参数 EBV 决定,该参数可以根据观察到的两条线的通量比计算得出,这两条线应该具有固定比率 r:

@pymc.deterministic
def EBV(Ha=Ha, Hb=Hb, r=r, R_V=R_V, Ha_l=Ha_l, Hb_l=Hb_l):
    kHb = gas_meas.calzetti_k(lams=np.array([Ha_l]), Rv=R_V)
    kHa = gas_meas.calzetti_k(lams=np.array([Hb_l]), Rv=R_V)
    return 2.5 / (kHb - kHa) * np.log10((Ha/Hb) / r)

我对 R_V 的值也有先验。

测量值本身表示为正态分布,例如

NII_6584 = pymc.Normal(
    'NII_6584', mu=f_row['[NII]6584'],
    tau=1./e_row['[NII]6584']**2.,
    observed=True, value=f_row['[NII]6584'])

我想获得 R_VEBVqz 的估计值。但是,当我从所有这些中创建 pymc Model 时,我被告知 Deterministic 对象不能具有观察值:

TypeError: __init__() got an unexpected keyword argument 'value'

首先,我是否误解了 Deterministic 对象的性质?如果是这样,我还能如何根据未直接观察到的值进行推断?

其次,我是否正确构建了观察结果?我必须将观察到的通量指定为均值和值参数,这似乎很奇怪,但我不清楚除了对通量均值和方差进行建模之外还能做什么,这似乎不必要地复杂。

如有任何建议,我们将不胜感激!

我认为您没有正确构建您的观察结果。这不是一个最低限度的工作示例,但也许我们可以消除一些困惑。

首先,我不认为 @deterministic 装饰器接受参数 value = <something>。不清楚您的哪些确定性陈述是实际模型,但请尝试将您的代码转换为以下模板:

#Define your randomly-distributed variables (I'm assuming they're normal)
q = pymc.Normal(name,mu=mu,tau=tau)
z = pymc.Normal(name2,mu=mu2,tau=tau2)

#Define how you think they generate your data
@pymc.deterministic
def NII_SII_th(q=q, z=z):
    return NII_SII_g(np.array([q, z])) #this fcn is defined somewhere else

#Your data array
f_row['[Nii]6584']=[...]

#Now link your model and your data
obs = pymc.Normal(modelname,mu=NII_SII_th,
                  observed=True, value=f_row['[NII]6584'])