Calculating highest density intervals requires multiple calls to calibrator

[joeroe]

Because we only cache the probability distribution, TAQ, and TPQ on cal records for c14 calibrations. I don't want to add HDI to that, because it belongs to specifically to (posterior) radiocarbon distributions, not to probabilistic dates more generally. Also, the choice of threshold is arbitrary and it would be nice to allow the user to select different ones (pending #328).

One solution would be to calculate the HDI in Ruby, perhaps building on https://github.com/joeroe/ruby-radiocarbon/, which would also make some related issues with calibrator go away (#323, #327, https://github.com/xronos-ch/xronos.config/issues/6).

On the other hand, the only actual problem with this—apart from the inelegance—is the possible performance hit, and I'm not sure that's significant. This is Cpp code we're talking about.

---

Incidentally, I went down a rabbit hole of terminology and algorithms in looking for a solution for this for ruby-radiocarbon. These were particularly useful reading:

• https://easystats.github.io/bayestestR/articles/credible_interval.html
• https://stats.stackexchange.com/questions/381520/how-can-i-estimate-the-highest-posterior-density-interval-from-a-set-of-x-y-valu (both answers)

As I understand it, what archaeologists refer to as "sigma ranges" are more properly the highest density intervals, regions, or sets.

Sometimes a p for posterior sneaks in there, but I think that's not quite right because they can be extracted from any kind of probability density distribution; p for probability makes more sense, but seems redundant. Some say that an 'interval' should be reserved for a single continuous range, but that seems a bit arbitrary to me, and 'interval' is more intuitive when we're talking about time.

[joeroe]

Actually I can think of some benefit of doing this outside calibrator:

• We could apply it to things other than radiocarbon dates. Although would assume we have another type of date described by a non-parametric function for which the HDI would be useful, so it's pretty hypothetical.
• We could experiment with different measures. Single HDI vs. multiple HDIs, for example. Or the easystats docs suggest that the 'shortest probability interval' is better (https://easystats.github.io/bayestestR/reference/spi.html)