Explicitly set reference epoch, save chi2 to object, and fixes plotting for stationary tracks

[wbalmer]

This is my attempt to both address issues #39 and #40 , but also to fix the stationary track component of trackplot.

This branch adds/updates:

• stationary_params, stationary_loglike, and stationary_chi2_red
• creates median_loglike and median_chi2_red once fit completes, or when load_results is called
• the ref_epoch_idx kwarg (defaults to 0), which sets ref_epoch explicitly during init
• corrects (I think) the ra0 and dec0 calculation so that the uniform prior boundary can be decreased or increased without affecting functionality
• updates stationtrackplot to reflect changes to trackplot
• removes ref_epoch in trackplot and takes directly from system.ref_epoch

I'm not in a hurry to merge this, but the stationary track issues have been bugging me for a while.

I began by thinking about the figure in Nielsen+2017 where they anchor the stationary track to the 2018 observation, and believe we should allow the user to set explicitly before the fit the observation from which to generate the stationary case.

Now, the user should be able to explicitly generate the stationary track plot starting from any given observation and determine the approximate chi2_r for that case, which can be useful when there are fewer observations. The forward and back propagation to determine the candidate position in absolute coordinates assumes the host star values from Gaia (e.g. it doesn't generate the track from the exact measurement, but from the projection of that measurement onto absolute coordinates based on the Gaia trajectory at the epoch of observation), so there is usually a small systematic between the start of the track and the actual measurement in the final plot, but I think this is okay at the level of accuracy someone would be interested in the stationary case (and of course this offset gets fit for when evaluating the non-stationary, finite case).

I'd appreciate @gotten 's review (and @tomasstolker if you have a chance!), and whether you think that setting this ref_epoch explicitly makes sense. Especially with regards to the chi2_red objects (for some reason, I'm not sure if I've calculated this correctly) and the ra0 dec0 calculations with SkyCoord.

[gotten]

Thanks for the PR. I'll take a more detailed look at it later this week.

Minimizing the distance to the observed offset and radecdists at a specific observed Epoch with RA, DEC at 2016 as the free parameters should do the trick as well. The free parameters would optimize to the ra0 and dec0 values needed for the stationary track via radecdists. The advantage would be that it takes into account the same assumptions and subtleties as the rest of the code (gravitational lensing by solar system bodies etc.). Maybe/hopefully the outcome will look the same. I have the code to do this somewhere.

[wbalmer]

I threw together a quick function in commit 111d685 that does that minimization of radecdist and the observation at the ref_epoch with scipy.optimize.minimize to determine ra0, dec0. It gets called in the place within system.init that ra0 and dec0 were originally defined, but the function itself might need to be cleaned up.

Fixing that small difference means the plots look a lot better!

[gotten]

The stationary track looks good, goes straight through the selected datapoint. We can still debate on the best way to set the start and end points of the tracks. The reference epoch can now be easily adjusted to an observed datapoint but because the start and endpoint depend on it (they depended on the midpoint/custom reference point before) they change with a new choice of reference datapoint. It works but might not be ideal when you are happy with the length of the tracks but just want the track to be pinned to a different datapoint. It's a minor thing that we can think about.

The chi2 in the PR is now log(-2*loglike). If the log likelihood is ~ -chi2/2 then the chi2=-2*loglike without the log ? I'm doing a run now to see how different the values are.

Do parallax, pm and dec of the BG star count as free parameters for the stationary track ? If not, we should use 8 for the number of free parameters instead of 11. So in pseudocode: -2*LL / (n_observations-n_dims+3) for the stationary track and -2*LL / (n_observations-n_dims) for the regular fit. n_observations has to include both RA and DEC data so instead of len(self.epochs) it will be 2*len(self.epochs) ?.

Edit: chi2_red=20.41 for stationary and 7.42 for median track (using 11 as number of parameters in both cases). Is the median track close to the maximum likelihood track ? Maybe the errors are underestimated or maybe the two other stars in the system contribute.

Edit2: Because we force the stationary track through 1 epoch that one doesnt contribute to the chi2. Maybe that means +1 instead of +3 in the denominator because we effectively lose 1 epoch with 1 ra and 1 dec measurement. Which slightly inflates that reduced chi2.

[wbalmer]

Right, thanks for catching that! I was coding without coffee when I wrote up that chi2 line... commit 465c0f3 has

self.stationary_chi2_red = -2.*self.stationary_loglike/((2*(len(self.epochs)-1))-self.ndim) and self.stationary_chi2_red = -2.*self.stationary_loglike/((2*(len(self.epochs)-1))-self.ndim)

does that look right?

[gotten]

OK, that seems to work.

I think we should give dummy_loglike a more descriptive function name because it seems to be far evolved from a loglike :).

I should have done the start and end point differently than a reference date and days_backward and days_forward but it evolved that way. Ideally we should move away from it because it was not very intuitive to begin with and in the current PR it is not a stable value because it depends on the choice of reference datapoint (although the default point is the same if not changed). We should think about what is the best / user friendliest way to implement similar functionality. Default behaviour: take the min and max time of the datapoints and add a certain fixed margin before and after ? With an option to override it based on decimal year (for instance setting [2010.0,2030.0] to some keyword parameter ? We can then drop the reference date and days_forward & backward in the associated code. The default scenario would ensure that if somebody just ran the code with minimal syntax they would still get a decent looking plot after waiting for 10+ minutes. Orbitize! uses start_mjd (in mjd) and sep_pa_end_year (in decimal years) parameters. Not ideal in terms of consistency. Their defaults translate to 2000 and 2025 at the moment. I'm open to suggestions, also from @tomasstolker because he has tweaked the days_forward/backward code a bit in past commits.

I haven't yet checked the code that produces the initial estimate from which to start the minimization but also acknowledge it doesn't have to be perfect to reach the minimum. The old version was inexact for multiple reasons. The current version is a big improvement as it goes to the required end result in two phases.

At which moment should we prune the code to remove unnecessary things ? Is that done before a 1.0.0 release or earlier in the development? We have some unused code snippets and maybe I tried to hard to have the code still work without many priors which clutters the code with if statements.

Because I'm not too experienced with development via github I wonder what is the best approach. Accept the PR and fix things later with new issues or try to improve things before pulling?

[wbalmer]

If the chi2 calculation looks good, and the stationary bug is somewhat alleviated, I'll merge this so we can keep things tidy and open a new issue to fix the implementation and defaults of the plotting timescales once and for all :) .

[gotten]

Sounds good. Thank you :)

[wbalmer]

Closes #39 and #40