joshspeagle / brutus

Modeling stellar photometry with "brute force" methods
MIT License
32 stars 5 forks source link

Including parallax during optimization #6

Closed joshspeagle closed 5 years ago

joshspeagle commented 5 years ago

Allow for the parallax to be included as part of the optimization step when solving for (s, Av).

Currently, the procedure I use to incorporate the parallaxes into the fit occurs in three steps, all after solving for (s_mle, Av_mle) and the respective covariances. In the first step, I pre-select models with good s_mle that match the provided parallax. In the second, I then integrate over the respective PDFs to evaluate the overlap integral using Monte Carlo draws (and importance reweighting) to get p_model. In the third, I draw samples by first drawing the model model | p_model and then the relevant samples from it s, av | model, weights.

Step 1 currently can remove a lot more models than needed, and Step 3 implicitly relies on good coverage, both of which fail if the photometry isn't great. If I can increase the coverage during the optimization step and then "subtract" its contribution using importance reweighting, that'd be a huge win-win overall.

joshspeagle commented 5 years ago

This has now been partially dealt with in a973b62535c4a23022825d44e732cb579a4e02c6. The initial step now includes (approximate) scale-factor errors. Turns out the importance reweighting didn't work out quite as I expected, but I found that increasing the number of draws from 150->250 helped shore up those results a bit.

In general, the actual distance uncertainties for a particular model tend to be quite small, so incorporating the Gaia parallax directly into the derived scale-factor doesn't actually help very much when we have many photometric bands. In the case where we don't, we then end up with a lot more models having equivalent weights, and increasing the number of distances realizations per model then allows us to marginalize over them more effectively when sampling distances anyway. So this effectively deals with both cases.

While I could still try and incorporate the Gaia parallaxes into the derived scale, av, and cov_sa values in the future, the code already is designed to spit out samples and re-process those with the provided priors anyway, so it's a bit of a moot point. Users should be grabbing the samples anyway in most cases, and I'll be updating the demos/docs to make that clear when the code development settles down.