Modeling separations

[cmlamman]

Looking at literature, ie https://arxiv.org/abs/1901.06364 and https://arxiv.org/abs/1209.4087, I'm worried that a gamma distribution is not right for modeling separations. It appears that the separation distribution is very broad and it doesn't seem like gamma distributions are good for that (but I could be wrong).

I'm thinking a Weibull or Lognormal

[cmlamman]

It seems like the most common choice is to use a gaussian of the log of separations. I feel like this should be an easy switch from our current model, however I've gotten stuck on how to properly sample from a distribution like this.

Here is a short notebook summarizing my problem: separation_distribution.zip

[cmlamman]

Update: figured out how to use the lognormal distribution. working on putting it into model

[vnmanoharan]

scipy.stats.lognormal is notoriously difficult to figure out. See some hints here: https://nbviewer.jupyter.org/url/xweb.geos.ed.ac.uk/~jsteven5/blog/lognormal_distributions.ipynb and here https://stackoverflow.com/questions/28700694/log-normal-random-variables-with-scipy

[cmlamman]

@vnmanoharan Thanks! I think I've mostly figured it out, here is the current fit I've gotten

And it's slightly better when accounting for the variable cutoff at low separations that depends on contrasts:

Do these look acceptable? I have yet to test this model on sets of simualted data with different parameters used to make the initial population (to ensure this model holds for different datasets)

[vnmanoharan]

If you can put a self-contained version in a notebook (with the model included in the notebook) in a new branch, I'll take a look.

[cmlamman]

Thank you, I've added a self-contained notebook which simulated and fits data, simulate_and_model.ipynb in the branch model_testing.

Sometimes the sampler runs out of bounds, but it seems that most of the time this notebook will run successfully (another issue)

[vnmanoharan]

I was not able to get it to sample without a mass matrix error. The problem seems to be in the likelihood for the inverse contrast ratios:

cr_observed_inverse = pm.Normal('cr_observed', mu=contrast_ratios_inverted, sigma=cr_err_inverse, observed=cr_inverse)

If I comment that out and the truncated_crs line I can get it to sample the separation distribution. Don't know what's happening. Maybe check that your uncertainty calculation is correct for cr_err_inverse?

[cmlamman]

Thanks, I think that was the issue! It ran 5/5 times when I used cr_err_inverse = 1 / cr_err instead of = cr_err / cr I think maybe usual error propagation isn't appropriate here

[cmlamman]

The fit to our simulated data looks good (see tutorial end), however the width of the separation distribution is very sensitive to the choice of prior. I've tried different priors (Weibull, lognormal, gamma, uniform) and even set the mean of this distribution to be a function of the width of the input angular separations. However if I simulate the same data but with a width of 0.5 or 1.5 (current is 1.0), the fit is always much worse.

Realistically, I don't think this width will be very far from 1 so maybe it's not worth worrying about?

(attempts for separation priors visible in comments in simulate_and_model.ipynb)

[vnmanoharan]

The actual uncertainties (assuming Gaussian on cr) should be cr_err/(cr**2). This also doesn't sample. I think the problem is that some of the cr are very small, and so taking the inverse leads to huge uncertainties for some data points. It may be more numerically stable to do the following:

• model the inverse mass ratios as Pareto distributed (you are already doing this)
• then define your Deterministic to go from inverse mass ratios to (non-inverted) contrast ratios
• then make your observed variable the (non-inverted) contrast ratios

In other words, the observations are the actual contrast ratios, but the underlying model is in terms of the inverse mass ratios. This way you don't have to convert uncertainties or convert data.

[cmlamman]

This was the first method I tried, but with no success, and going back to it I still can't get it to run without sampling errors (even just fitting for contrasts and no truncation at all).

Maybe the issue is that the power index cannot be <1, because that doesn't match the contrast ratio distribution at all (it goes to infinity at 0). However it can be <1 and still match the inverse contrast ratios well. So this method requires a steep cutoff around 1, but that doesn't work with the inverse contrast ratios.

Maybe I have to use inverse contrast ratios the entire way through or not at all

Update: I was able to get it to run consistently using inverse contrast ratios the entire way through and the proper error formula. Perhaps the issue before was not that the errors were to big, but that I had a prior which allowed the power law to drop below 1

[cmlamman]

Turns out its not running as consistently as I thought, this is still an issue. To summarize this problem, the three methods we can do all seem to have issues:

1. Fit to a power law of mass ratios: problem is that you have a sharp cutoff in the distribution at 1
2. Fit to a Parento distribution of inverse mass ratios but then transform to contrast ratios for likelihood: I think the issue here is that you need to put a hard prior so power_index > 1 or else the sampler runs out of bounds. But if you do this then the hard cutoff causes problems in the Parento distribution
3. Fit to a Perento distribution of inverse mass ratios and calculate likelihood with inverse contrast ratios: problem is that the errors for inverse contrast ratios are huge.

It seems like the simplest option for now would be to just scale / maybe take the sqrt of contrast ratio errors and go with option 3. This is the only way that runs all the time, and it produces a pretty good fit. I still don't like to improperly handle error propagation though.

[cmlamman]

Closed, am confident in how the model used lognormal now