pymc: Inferring parameters based on functions of observables

Question

I have observations of several optical emission lines, and I have a model that predicts several (flux) ratios of those lines, based on two parameters, q and z, which I want to infer.

I have created @pymc.deterministic objects that take values of q and z (each of which has uninformative priors over some physically-interesting region), and turn them into a "predicted" ratio. There are about 7 ratios, and they have the form:

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

I can also define the ratios derived from observations, such as

@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))

where, for instance, NII_6584 is the observed flux of one of the lines and rcf_NII_6584 is the flux correction for that same line. These corrections are themselves determined by the line wavelengths (known with infinite precision), and by a parameter EBV, which can be calculated from the observed flux ratio of two lines that are supposed to have a fixed ratio 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)

I also have a prior on the value of R_V.

The measurements themselves are expressed as Normal distributions, such as

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'])

I would like to get estimates of R_V, EBV, q, and z. However, when I make a pymc Model from all these, I am told that Deterministic objects cannot have observed values:

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

First, am I misunderstanding the nature of Deterministic objects? If so, how else do I infer based on values that are not directly observed?

Second, am I constructing the observations correctly? It seems odd that I'd have to specify the observed flux as both the mean and the value argument, but it's not clear to me what else to do, other than also model the flux means and variances, which seems unnecessarily complicated.

Any advice would be appreciated!

score 2 · Accepted Answer · answered Mar 29 '16 at 19:45

I don't think you're constructing your observations correctly. This is not a minimum working example, but maybe we can clear up some confusion.

First off, I don't think the @deterministic decorator takes an argument value = <something>. It's not clear which of your deterministic statements is the actual model, but try to translate your code into the following template:

#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'])

This is fairly helpful, so I'm accepting, since it's the only attempted answer I've gotten. Thanks for your help! — DathosPachy, Mar 30 '16 at 13:55

pymc: Inferring parameters based on functions of observables

1 Answers1