Fix NA density value when a bound coincides with an extremum of x

[haomeng797-ship-it]

Closes #269.

The cause

density_bounded() trims by interpolating the reflected density onto seq(min(x), max(x), length.out = n):

x_trimmed = seq.int(range_x[[1]], range_x[[2]], length.out = n)
f = approx(d$x, f, x_trimmed)$y

d$x is a slice of the wider grid seq(bounds[[1]] - width, bounds[[2]] + width, length.out = n_unbounded), so its endpoints reproduce bounds only up to floating point error. When a bound happens to equal the corresponding extremum of x, that error is enough to break the interpolation: x_trimmed ends exactly at max(x), d$x ends a couple of ulp below it, and approx() returns NA for that last point.

In @ASKurz's reprex the discovered upper bound is exactly max(x) = 1, and the grid lands at 0.99999999999999978, short by 2.22e-16. That single NA then reaches quantile() inside hdi(), which is where the error surfaces:

Error in quantile.default(dist_y, probs = 1 - .width, names = FALSE) :
  missing values and NaN's not allowed if 'na.rm' is FALSE

A smaller deterministic reproducer, no sampling involved:

density_bounded(c(0.5, 1, 1.5, 2), bounds = c(0.1, 2))$y |> anyNA()
#> [1] TRUE

The fix

approx(..., rule = 2). The trim branch only runs when bounds already contains range(x), so anything outside the interpolation range is floating point noise at the endpoints, and constant extrapolation over a few ulp is exact to the precision available. This leaves the returned x grid unchanged, so trim still means "the output spans exactly the range of the data".

I considered clamping x_trimmed to range(d$x) instead, but that would make the returned grid stop just short of max(x) and quietly weaken what trim promises.

Tests

Added a regression test using the deterministic reproducer above. It fails on main with both the NA and the hdi() error, and passes with the fix. test.density.R and test.point_interval.R pass; the remaining failures in the full suite are the pre-existing vdiffr snapshot differences from #272, which also fail on a clean checkout here.

One thing I noticed while testing, not addressed here since it is a separate question: on the reprex's p^3 posterior, whose true density is monotone, hdi() now returns two disjoint intervals rather than one. That looks like ordinary noise in the density estimate creating a shallow local minimum, not a bug in this code path, but it may be worth a look if you think the density should be smooth enough there to avoid splitting.