Handling Gaia DR3 errors for the host star

[gotten]

Considering the Gaia DR3 position, pm, parallax and rv for the host star have errors that become relevant at GRAVITY's accuracy we should implement a way to include the Gaia error into the fit. For instance, HD 131399A has a positional error of 30 microarcsec in RA and DEC. Otherwise the posteriors for the background object do not include the error in the host star.

A way to implement it, is to redraw the host_star location from the DR3 distribution. This is computationally intensive and will at least double the execution time because we will need to setup the catalog and do transform_cat for the host for every new location.

Not sure how to implement it efficiently in combination with a fitter like emcee / dynesty. Give the host star parameters as additional parameters to the model and add the likelihood function for Gaia DR3 ? This would double the amount of parameters that are being optimized which might lengthen the convergence time and further increase the amount of calls needed to the model ?

Here I can also say that for the purpose of checking if a point source follows a background-like trajectory this is not needed. But it is needed for accurate BG star parameters. Not an immediate priority.

[tomasstolker]

Doubling the parameters should still be fine I think, so would be interesting to test! Currently, the fit takes about 5 min for HD 131399 on my laptop with dynamic=False.

From what I quickly tested, both make_cat_entry and transform_cat take about 20 microsec so that shouldn't be a problem?

For the implementation, I think you can indeed just include these additional free parameters, and sample the host star parameters in prior_transform with transform_normal for a normal distribution.

[tomasstolker]

Depending on the number of parameters, there are different sampling methods in dynesty: https://dynesty.readthedocs.io/en/latest/quickstart.html#sampling-options

@wbalmer added the sample_method parameter to fit, which is the sample parameter in dynesty. You could try set it to 'rwalk' to test if it helps with the convergence when including the extra host star parameters.

[gotten]

A small search tells me that there is no PPF for a multivariate normal distribution, so we have to approximate it by N independent normal distributions without a covariance term. Not ideal but better than nothing.

PDF, CDF and likelihood does have a definition for the multivariate normal, so an emcee implementation might work?

Would it be possible to use uniform priors and use the loglikelihood to have the behavior of a multivariate prior ?

Edit (not sure if relevant): https://stats.stackexchange.com/questions/436817/sampling-prior-covariance-matrices-nested-sampling https://www.sciencedirect.com/science/article/pii/S0047259X06000509 https://johannesbuchner.github.io/UltraNest/priors.html#Dependent-priors most promising: https://pints.readthedocs.io/en/latest/log_priors.html#pints.MultivariateGaussianLogPrior https://github.com/pints-team/pints/blob/728a22310447f3768629358427a43492007777f5/pints/_log_priors.py#L819 https://github.com/pints-team/pints/blob/728a22310447f3768629358427a43492007777f5/pints/_log_priors.py#L984

Edit 2: I'm reasonably optimistic that the pseudo-inverse approach works. Tested it with 2D multivariate gauss with as input uniformly drawn numbers.

mean=[5,-3] cov=[[1, -0.5],[-0.5, 1]] 2D histogram:

[wbalmer]

These PPFs are tricky! The pseudo-inverse seems like a reasonable approach, I believe that scipy uses the pseudo-inverse to calculate the CDF for their multivariate_normal function.

[gotten]

A teaser with dummy covariance matrix. To be pushed to fork when it is closer to implemented In the class backtracks.System with ra,dec,pmra,pmdec,parallax,rv of the host star:

self.HostStarPriors = HostStarPriors(np.r_[self.rao,self.deco,self.pmrao,self.pmdeco,self.paro,self.radvelo],np.eye(6))

with the class HostStarPriors defined in utils.py and imported into backtracks.py In an interactive python session:

import backtracks import numpy as np b=backtracks.System("HD 131399A",candidate_file="./tests/scorpions1b_orbitizelike.csv") randomx=np.random.random([2,6]) # unit samples b.HostStarPriors.transform_normal_multivariate(randomx)

we get an array drawn from our fake prior

array([[223.16773393, -33.64350581, -31.72359809, -31.01404474,
         10.63760254,  14.33849705],
       [223.29501983, -33.08704849, -30.62811149, -30.28465357,
         10.81069543,  16.36160253]])

The code uses a class because the distribution doesn't change so part of the maths can be done in advance. The evaluation of the PPF can then do the rest of the maths quicker.

Now the rest of the owl !

Edit:

it converged !: iter: 2431 | +100 | bound: 124 | nc: 1 | ncall: 22280 | eff(%): 11.411 | loglstar: -inf < -100.245 < inf | logz: -123.701 +/- 0.475 | dlogz: 0.005 > 0.500

[gotten]

Plots ! I'm suppressing the dimensions of the host star by default.

I think the neighbourhood plot still crashes. Needs a small bugfix.

Do you think the plots make sense? No idea how to interpret the Evidence plot.

Edit: neighbourhood plot fixed.