sczesla / PyAstronomy

A collection of astronomy-related routines in Python
152 stars 35 forks source link

Make KeplerEllipseModel work with fitEMCEE method #25

Closed elnjensen closed 7 years ago

elnjensen commented 7 years ago

I found that the KeplerEllipseModel code would not work well with the fitEMCEE method, because it would try to evaluate orbits with non-physical orbital parameters, such as e < 0 or e > 1. Even if restrictions are set on these parameters with setRestriction, these restrictions are enforced only after the chi-squared evaluation (via the penalty factor), so the code would still try to evaluate the orbit first, returning an error, and fitEMCEE would fail, reporting ValueError: lnprob returned NaN.

The solution is twofold. First, implement limits on these parameters (e.g. eccentricity and various angles in the orbit) via defining priors. Following the recommended emcee method for parameter limits, these priors return -numpy.inf for the log probability if parameters are out of bounds. Second, I have modified the fitEMCEE code slightly so that it checks the priors before trying to evaluate a given orbit, and if the prior sum is -numpy.inf, it short-circuits and simply returns that value. In addition to avoiding unphysical parameters, this should also be more efficient overall, avoiding unnecessary function evaluations in cases where the prior probability is zero.

Unrelated to the above, I found that the code for fitEMCEE was setting the fitting object to the lowest probability solution rather than the highest at the end of a run, so I fixed that.

If you think the above changes are useful, I could also make some modest changes to the docs for fitEMCEE to show how to set parameter limits in this way.

sczesla commented 7 years ago

Hi,

Thanks for the update! Indeed, I think these changes are very useful.

I completely agree with your diagnosis of the problem with KeplerEllipse and fitEMCEE (and many other models), and your solution is great. Ideally, I think, there should be a kind of built-in prior in the model, taking care of this problem, but that may be a different story.

Incorporating this into the documentation will be useful. This appears to me a rather frequent problem, which will also pop up when other models with similar boundary conditions are used. Any suggestions will certainly be welcome!

Thanks for fixing the "lowest probability bug" (it truely amazed me).

Cheers, Stefan