pymc:根据可观察对象的函数推断参数
pymc: Inferring parameters based on functions of observables
我观察到几条光发射线,我有一个模型可以根据我想要的两个参数 q
和 z
预测这些线的几个(通量)比率推断。
我创建了 @pymc.deterministic
个对象,这些对象采用 q
和 z
的值(每个对象在一些物理上感兴趣的区域都没有信息先验),并将它们变成 "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_V
、EBV
、q
和 z
的估计值。但是,当我从所有这些中创建 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'])
我观察到几条光发射线,我有一个模型可以根据我想要的两个参数 q
和 z
预测这些线的几个(通量)比率推断。
我创建了 @pymc.deterministic
个对象,这些对象采用 q
和 z
的值(每个对象在一些物理上感兴趣的区域都没有信息先验),并将它们变成 "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_V
、EBV
、q
和 z
的估计值。但是,当我从所有这些中创建 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'])