GP refactor

[sblunt]

I made a reasonably large change to the GP functionality in RadVel. Here's the gist:

When I first wrote the GP code, I treated data from each instrument as a separate realization of a Gaussian Process. This meant that although the GPs could share hyperparameters (e.g. they had the same stellar rotation period, etc), I assumed that one data point taken by HIRES followed immediately by one data point taken by HARPS were uncorrelated. This was something that Molly Kosariak and Trevor David pointed out as being a not great assumption! I've been meaning to fix this for a while.

The new implementation assumes all RVs are correlated with one another, although the GP amplitudes are still allowed to be different. The GP prediction for a given instrument is now just a multiple of the GP prediction for another instrument. This ended up being a significant refactor under the hood, which primarily involved updating the GPLikelihood object to inherit from the CompositeLikelihood object rather than the RVLikelihood object. The only change command-line users will see is that now hnames, defined in the config file, is a list instead of a dictionary. I've made this change in the example GP config file, and it runs fine. I wrote this config file to reproduce Fei's fit (which, in his 2017 paper, assumes the data have the same GP amplitude), and it still does very well.

So, are all past RadVel fits using GPs wrong? No, but they are possibly less precise than they would be using this new assumption (i.e. posteriors might be a little wider). I essentially just left out some information in the data when computing the GP posteriors. In addition, if you did a GP fit using data taken by only one instrument, this branch will perform identically. I think it would be good to clarify this in a paper that folks can reference-- I'm currently writing up a GP fit, and I think it would be good to include an appendix writing out this point with all the gory math.

Still to do: I've only updated the QuasiPeriodic kernel. I saw @vandalt opened a PR adding some support for other celerite kernels, and was wondering if you would be willing to rebase on this branch (as you mentioned in your PR), since I think that would save us both a lot of extra work. Let me know what you think!

[vandalt]

Sounds good! I'll mark #358 as draft and rebase it on this branch.

Once this PR gets merged, will it still be possible to use separate GP likelihoods and then combine them in a CompositeLikelihood (or even in a "GPCompositeLikelihood") ? While this is indeed not the best approach for multiple RV datasets, I feel like it could still be valuable for custom model (for example a separate GP with shared parameters for a light curve) and/or for instrumental effects specific to one dataset.

[sblunt]

Awesome, thanks!! No, under this refactor that wouldn't be possible. We certainly could add it back in though. For joint photometry training, we can use the NumericalPrior to set posteriors from a GP trained on photometry as priors for the RV GP, so I don't think we'd need it for that. However, the individual instrumental effects modeling is a good idea. It would be good to add that capability back in on top of the current implementation. Since that's still not really supported in the current master version of the code, though, I'd prefer to focus on finishing up this implementation first. Does that make sense to you?

[sblunt]

Hm, thinking more about it, I think I was falling into the same trap with my thinking about the joint fitting of photometry and RVs vs training on photometry then using those posteriors as priors. These are only equivalent if we assume the photometry and RVs are uncorrelated, which isn't true when we're doing GP fits.

[sblunt]

If we'd want to do a true joint fit, I think we'd want to use something like Rajpaul's FF' method, which tries to explicitly model the relationship between photometry and RVs. Otherwise, we might be overspecifying the nature of the correlation between them.

[vandalt]

It would be good to add that capability back in on top of the current implementation. Since that's still not really supported in the current master version of the code, though, I'd prefer to focus on finishing up this implementation first. Does that make sense to you?

Yes! I could try to implement this when the current PR is merged.

For the photometry or activity indicator modelling, I think modular GPs that can be combined would be useful for a simple joint fit or for something more involved like the framework from Rajpaul+ 2015. Users would then have the ability to test different methods, even though the Numerical prior option is likely good enough for most scenarios.

[sblunt]

Sweet, sounds great!

[bjfultn]

Should I merge this @sblunt ?

[sblunt]

Not yet @bjfultn, I need to refactor the other kernels.