Assertion failed in _linear_interpol

[mentics]

The implementation of _linear_interpol seems to be numerically unstable. Using Float64's, I ran into this case:

x_lo = 0.049999999999999996
x_hi = 0.050199999999999995
y_lo = 0.9729537570597172
y_hi = 0.9733565551446345
x = 0.05
EmpiricalDistributions._linear_interpol(x_lo, x_hi, y_lo, y_hi, x)

This results in the @assert throwing because y = 0.9729537570597171 which is just under (last digit) y_lo.
I'm not sure of the optimal solution (and maybe different cases would choose different tradeoffs). I found some discussions of it:

https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2018/p0811r2.html
https://math.stackexchange.com/questions/907327/accurate-floating-point-linear-interpolation

and many more.
In this case, I think you only need to handle x in [x_lo, x_hi], so it can be a simpler than the lerp in that first link.

I don't know what the definitive solution is, but I used this and it seemed to work for my case:

function _linear_interpol(x_lo::Real, x_hi::Real, y_lo::Real, y_hi::Real, x::Real)
    if x == x_hi
        return y_hi
    end
    T = promote_type(typeof(y_lo), typeof(y_hi), typeof(x))
    w_hi = T(_ratio(x - x_lo, x_hi - x_lo))
    y = y_lo + w_hi * (y_hi - y_lo)
    @assert y_lo <= y <= y_hi
    y
end

It's not perfect. I think it's possible for y > y_hi, but it hasn't hit the assert for me yet (running it on millions of data points), and if it did, I think I'd just clamp it (if y > y_hi, y = y_hi)

[oschulz]

Hm, how accurate do we need to be?

[mentics]

Right, some of the answers in that second link went much farther than is typically needed, and I wouldn't want to pay the cost of that level of accuracy. To me, the important thing is satisfying the specific properties that will avoid potential issues with downstream calculations. And not throwing incorrect exceptions, of course :)

The first link has an implementation that seems reasonable, except I think it supporting extrapolation which we don't need here. I'm not an expert, but just something to "tighten up" the ends of what you have is probably sufficient.

[oschulz]

Want to prepare a pull request?