Modeling the ratio of two stochastics
Todd
I'm wondering if someone has some suggestions for how to model the
ratio of two stochastics with PyMC. Specifically, I have observations
(with associated Poisson counting errors) of the flux (counts/sec)
from quasars in two ultraviolet wavebands, and I'd like to fit a
lognormal function to the distribution of the flux ratios, r = f1/f2.
I have a simple model for inferring the fluxes in each band, assuming
that all the source and background counts are well-described by
Poisson distributions:
f_obj = pymc.Uniform('f_obj', 0., 100., value=tot_data - sky_data)
f_sky = pymc.Uniform('f_sky', 0., 100., value=sky_data)
@pymc.deterministic
def f_tot(f_obj=f_obj, f_sky=f_sky):
return f_obj + f_sky
d_sky = pymc.Poisson('d_sky', f_sky, value=sky_data, observed=True)
d_tot = pymc.Poisson('d_tot', f_tot, value=tot_data, observed=True)
M = pymc.MCMC([f_obj, f_sky, f_tot, d_sky, d_tot])
This simple model works very well for estimating the observed fluxes,
including the uncertainty in the background. I'd like to feed the
posterior distributions for f_obj for the two bands for all the
quasars back into PyMC to infer the parameters of the lognormal
distribution, but I'm not sure how best to do it. I know that I can
construct the posterior distribution for the flux ratio for each
quasar by either directly using the posterior samples of the two
fluxes or by marginalizing the joint distribution of the fluxes:
p(r|...) = \int f2 p(f2*r, f2)df2
With the posterior distribution for the flux ratio, I suppose I could
write my own pymc.stochastic function for the observed ratio r and use
a model something like the following:
r_true = pymc.Lognormal('ratio', mu, tau)
r = myStochastic(...)
But I'm wondering if there's some more elegant solution that could
infer the ratio and the fluxes simultaneously by adding elements like
the following to my model for inferring the fluxes:
@pymc.deterministic
def ratio(f1=f1, f2=f2):
return f1/f2
{ratio} = pymc.Lognormal('varName', mu, tau)
where {ratio} is the set of flux ratios for all the quasars in the
data set. The problem I run into naively trying to pursue such a
solution is that I want "ratio" to be both deterministic (i.e.,
determined by the fluxes) and an observed stochastic where the
"observed" values change with each iteration of the sampler. Could I
write my own stochastic function with a "random" method that updates
the value of the stochastic at each iteration with the deterministic
ratio (instead of some truly random value)?
I'd be grateful for any advice.
Thanks,
Todd
Anand Patil
On Mon, Jul 13, 2009 at 6:30 AM, Todd<todd_...@mac.com> wrote:
>
> This simple model works very well for estimating the observed fluxes,
> including the uncertainty in the background. I'd like to feed the
> posterior distributions for f_obj for the two bands for all the
> quasars back into PyMC to infer the parameters of the lognormal
> distribution, but I'm not sure how best to do it. I know that I can
> construct the posterior distribution for the flux ratio for each
> quasar by either directly using the posterior samples of the two
> fluxes or by marginalizing the joint distribution of the fluxes:
>
> p(r|...) = \int f2 p(f2*r, f2)df2
The reason this is confusing is the requirement that the ratio have a
lognormal posterior is inconsistent with your model. The 'right' thing
to do would be to make a histogram of the posterior of the ratio using
the MCMC samples, and present that in place of the lognormal
distribution.
You can still present a lognormal approximation to the posterior, but
it's just going to be an approximation, so you may as well keep it
simple. You could estimate the parameters by simply taking the mean
and variance of samples of log(r), and then compute the
Kullback-Liebler divergence or Kolmogorov-Smirnov test statistic or
something to see how good of an approximation it is.
If you really want to force the posterior to be lognormal 'inside' the
model, MCMC is probably not the right tool for the job. Consider a
variational approximation.
Finally, I don't know how these things work, but it's worth having a
careful think about whether the f's for the two bands should really be
modeled independently. For example, should the two f_sky's be
correlated? If so, it could have a noticeable effect on the posterior
of the ratio.
Cheers
Anand
Anand Patil
Wait- do you want to summarize the posterior of a single ratio as
lognormal, or model the distribution of many ratios as lognormal? If
the latter, scratch everything I said. :)
The easiest thing to do is make ratio a lower-level variable than f1
or f2, and put both bands of all your datapoints in the same model
mu_ratio = Normal('mu_ratio',...)
tau_ratio = Gamma('tau_ratio',...)
ratios = Lognormal('ratio',mu_ratio,tau_ratio)
f1 = Uniform('f1',0,100,value = ones(n))
f2 = f1 * ratios
...
and then you get posteriors for mu_ratio and tau_ratio, as you (may
have) wanted.
Anand
Todd
On Jul 13, 6:36 am, Anand Patil <anand.prabhakar.pa...@gmail.com>
wrote:
>
>
>
>
>
> Patil<anand.prabhakar.pa...@gmail.com> wrote:
> > Todd,
>
a lognormal. I need to reconfirm this, but I think I quickly tried a
model like the one you suggest over the weekend and had trouble with
the induced correlation between "ratios" and f1. mu_ratio and
tau_ratio would actually be estimated correctly, but f1 would end up
with way too narrow a posterior distribution (i.e., much smaller than
the square root of the counts). Now that you've suggested a similar
model, I have more confidence that I'm on the right track. I'll try
again, perhaps with an AdaptiveMetropolis step method, and report back
on how things work.
Thanks again for your help!
Todd
Anand Patil
with unkown mean vector and covariance matrix. The predictive
distribution of f1/f2, conditional on those parameters, would then be
lognormal.
Todd
On Jul 13, 12:42 pm, Anand Patil <anand.prabhakar.pa...@gmail.com>
wrote:
>
> > On Jul 13, 6:36 am, Anand Patil <anand.prabhakar.pa...@gmail.com>
> > wrote:
> >> On Mon, Jul 13, 2009 at 10:51 AM, Anand
>
> >> Patil<anand.prabhakar.pa...@gmail.com> wrote:
> >> > Todd,
>
> >> > On Mon, Jul 13, 2009 at 6:30 AM, Todd<todd_sm...@mac.com> wrote:
>
> >> >> This simple model works very well for estimating the observed fluxes,
> >> >> including the uncertainty in the background. I'd like to feed the
> >> >> posterior distributions for f_obj for the two bands for all the
> >> >> quasars back into PyMC to infer the parameters of the lognormal
> >> >> distribution, but I'm not sure how best to do it. I know that I can
> >> >> construct the posterior distribution for the flux ratio for each
> >> >> quasar by either directly using the posterior samples of the two
> >> >> fluxes or by marginalizing the joint distribution of the fluxes:
>
> >> >> p(r|...) = \int f2 p(f2*r, f2)df2
>
> >> lognormal, or model the distribution of many ratios as lognormal? If
> >> the latter, scratch everything I said. :)
>
> >> The easiest thing to do is make ratio a lower-level variable than f1
> >> or f2, and put both bands of all your datapoints in the same model
>
> >> mu_ratio = Normal('mu_ratio',...)
> >> tau_ratio = Gamma('tau_ratio',...)
>
> >> ratios = Lognormal('ratio',mu_ratio,tau_ratio)
> >> f1 = Uniform('f1',0,100,value = ones(n))
> >> f2 = f1 * ratios
>
> >> ...
>
> >> and then you get posteriors for mu_ratio and tau_ratio, as you (may
> >> have) wanted.
>
> >> Anand
>
>
> > a lognormal. I need to reconfirm this, but I think I quickly tried a
> > model like the one you suggest over the weekend and had trouble with
> > the induced correlation between "ratios" and f1. mu_ratio and
> > tau_ratio would actually be estimated correctly, but f1 would end up
> > with way too narrow a posterior distribution (i.e., much smaller than
> > the square root of the counts). Now that you've suggested a similar
> > model, I have more confidence that I'm on the right track. I'll try
> > again, perhaps with an AdaptiveMetropolis step method, and report back
> > on how things work.
>
> > Thanks again for your help!
>
> > Todd
Your suggested model works beautifully!
Thanks!
Todd