More flexible error distributions

[Parrazyte]

As of now, the library doesn't allow non-symmetric uncertainties, or values at a boundary (in which the uncertainty is one sided). These cases are extremely common for e.g. astrophysics fitting results. I've done it myself locally and it's a very minor change. I think a proper implementation up to the standards of such an important library would be very beneficiary to the community.

[privong]

Thanks for the request. Can you please provide some more information on what you'd specifically like to see implemented? As you've likely seen, the code currently assumes Gaussian uncertainties so that it can draw new values. Implementing the changes your request would require generalizing that code, but making it truly general is probably not trivial. I could use more information on what you're looking for.

• non-symmetric uncertainties: Can you say more about what situations you mean? There are two situations that come to mind, one where the uncertainties are only apparently (but not actually) asymmetric and one where they are

  • uncertainties in log10(value) Quantities that have symmetric uncertainties will have asymmetric uncertainties in log10(). New values could be drawn by transforming the quantities back into their non-log values, drawing from a Gaussian, and taking log10() of the draw/perturbed value.
  • arbitrary posterior probability distributions: If the uncertainties are asymmetric because the posterior distribution is non-Gaussian, correctly drawing perturbed values requires passing the full posterior probability distribution and using that.
  • other: Can you provide an example of a parametric uncertainty distribution that you've encountered? In principle, the code could be generalized to handle uncertainties that are skewed/asymmetric Gaussians, but in those cases it would be best to explicitly pass values for the sigma, skew, kurtosis, etc..

• values at a boundary: Just to confirm, here you don't mean upper or lower limits, correct? Can you provide an example of a situation where this is the case?

[Parrazyte]

For your first point, I was referring to the second case of non-symmetric posterior distributions. I agree that parsing the whole distribution is better, but in e.g. my case, computing the posterior distribution is computationally expensive and adding this input would require significant changes in your code. I was then proposing to add the "ad-hoc" choice of splitting the distribution on each side (e.g. providing 2 values for the sigma which would then be passed for different widths of each half-gaussian).

For the values at a boundary, yes, indeed, this amounts to upper limits. Direct example since you're also doing some x-ray astronomy : Fitting results with xspec may peg parameters to their respective max/min values, in which case the uncertainties and posterior distribution from these parameters will only have a tail on one side (of course assuming that this tail is gaussian is very crude).

[privong]

Thanks for the quick reply. Regarding non-symmetric distributions, both implementing the draw from the posteriors and hacking together a split distribution are doable (and would require additional bookkeeping). I'd like to see or do some tests to see how robustly the latter approach recovers the "right" answer, comparing to the full/general distributions before implementing something like this though.

For the boundary values, as you describe them, I think that the correct way to handle them is to flag them as upper or lower limits. Perturbing them around in the one-side of the tail would not fully capture the fact that the model fitting could not constraint them being beyond the boundary (at least in the case you noted). If the boundary is a physical one (e.g., causality constraints, as in some pulsar system timing), then treating them as fully constrained values makes sense. But then it'd be preferable to use the full posterior distribution as above.

[privong]

As a follow-up, if it's too computationally expensive to get the full posterior distributions, how are you estimating the asymmetric uncertainties that you'd like to feed into pymccorrelation? Implementing an ad hoc method would probably benefit from trying to "match" the method of estimating the uncertainties when doing the point perturbation.

[Parrazyte]

Thanks for your answer. Your make a good point for the need for "physicality" for upper limits. Currently I estimate the asymmetry from errors computations at a given sigma on these parameters.

[privong]

Currently I estimate the asymmetry from errors computations at a given sigma on these parameters.

Can you say a bit more about this? How are you deriving/obtaining the sigma values?

[privong]

@Parrazyte Pinging you about my recent question. Can you please explain how you estimate the asymmetric uncertainties and provide a minimum standalone example? This would help me see how to implement such a capability into this library.

[privong]

@Parrazyte At this point I don't have enough information to address your requested changes. Given it's been a few months without more info, I'm going to close the ticket. If you'd like to provide more information, please feel free to repoen it.