exoplanet-dev / exoplanet

Fast & scalable MCMC for all your exoplanet needs!
https://docs.exoplanet.codes
MIT License
206 stars 52 forks source link

Consider whether exoplanet should use DensityDist rather than Potential #70

Closed mrtommyb closed 4 years ago

mrtommyb commented 4 years ago

PyMC3 recommended way to include a custom likelihood function is with pm.DensityDist, e.g.

loglike = pm.DensityDist('loglike', logp, observed=y)

However in the exoplanet examples, the the GP likelihood is implemented as pm.Potential. The results are the same but it might be sensible to implement the samples using the standard PyMC3 way of doing it. Advantages include clear separating out of the observables. Additionally, some of the nice plotting functions like model_to_graphviz ignore potentials.

Here's an example of what it would look like

pm.DensityDist("loglike", gp.log_likelihood, observed=(y - light_curve - mean))
dfm commented 4 years ago

I think this would be cool, but it would probably be better to separate y from mean + light_curve since that is actually the mean. I've been working on a better GP interface, but who knows how long it will take so I'd love a pull request to dfm/exoplanet-docs that changes the notation.

mrtommyb commented 4 years ago

Ok. I will try to work on it.

dfm commented 4 years ago

@mrtommyb: I implemented this on the flight to AAS. It's on the new_gp branch if you want to take a look: https://github.com/dfm/exoplanet/blob/new_gp/src/exoplanet/gp/celerite.py#L117

The syntax is now something like:

with pm.Model():
    kernel = xo.gp.terms....
    gp = xo.gp.GP(kernel, x, yerr**2, ...)
    gp.marginal("likelihood", observed=y)
bmorris3 commented 3 years ago

Hi @dfm,

I'm working on a model which uses exoplanet to compute a transit+phase curve+eclipse model. I have a derived parameter (an eclipse depth) produced in my phase curve model which is wrapped in pm.Deterministic, and I would like to enforce a prior on this variable (eclipse depth), in order to be consistent with a value in the literature. I was thinking the easiest way to do this would be to run something like:

        pm.Potential('spitzer_depth_potential', pm.Normal.dist(0, 300).logp(spitzer_depth))

where spitzer_depth is the Deterministic, and its measured value is 0 ± 300 [ppm]. Will this work if I'm running gp.marginal later (like you eventually settled on) to create define the likelihood? Or is this a bad idea?

dfm commented 3 years ago

@bmorris3: There shouldn't be any issue with that syntax, but a simpler one might be:

pm.Normal("spitzer_depth_prior", mu=spitzer_depth, sigma=300, observed=0.0)

Which should work equivalently and I find somewhat more readable. (You can think of this "prior" as, instead, "data" from a previous experiment.)

bmorris3 commented 3 years ago

Thanks for this, that's definitely more readable!