JuliaStats / Distributions.jl

A Julia package for probability distributions and associated functions.
Other
1.1k stars 414 forks source link

`0 <= cdf(truncated(Normal(2.5, 0.2), lower=0.0),x)` not always true. #1854

Closed lrnv closed 4 months ago

lrnv commented 5 months ago

Followup to #1853

Consider the following:

using Distributions
D = truncated(Normal(2.5, 0.2), lower=0.0)
cdf(D, 3.741058503233821e-17) > 0  # false
cdf(D, 1.4354474178676617e-18) > 0 # false
cdf(D, 8.834854780587132e-18) > 0  # false

I think I am simply hitting catastrophic cancellation in this line.

I understand that such numerical incauracies are always possible at borders, but would that be considered a bug ?

devmotion commented 5 months ago

I'd consider a cdf that returns values < 0 or > 1 a bug. This is the case in your example:

julia> cdf(D, 3.741058503233821e-17)
-3.00686029885372e-50

julia> cdf(D, 1.4354474178676617e-18)
-3.00686029885372e-50

julia> cdf(D, 8.834854780587132e-18)
-3.00686029885372e-50

I'm leaning towards clamping the result in cdf because my initial feeling is that the truncated "wrapper" can be expected to be based on the cdf of the untruncated distribution.

Alternatively, one could define cdf(d::Truncated, x::Real) = exp(logcdf(d, x)). This would also fix the example:

julia> exp(logcdf(D, 3.741058503233821e-17))
0.0

julia> exp(logcdf(D, 1.4354474178676617e-18))
0.0

julia> exp(logcdf(D, 8.834854780587132e-18))
0.0

But I'm not sure if it's a good default since in many cases this will be less efficient - even though the current result from cdf is accurate enough.

lrnv commented 5 months ago

I think clamping on the line that I quoted is better than cdf(d::Truncated, x::Real) = exp(logcdf(d, x)) for performance reasons. The same thing should be done to ccdf(d::Truncated,...) for the same reason

burtonjosh commented 4 months ago

Just came across this issue myself with the following values:

julia> cdf(
         truncated(Normal(2.122039143928797, 0.07327367710864985),
           lower = 1.9521656132878236,
           upper = 2.8274429454898398,
         ),
         2.82,
       )
1.0000000000000002

Julia version 1.10.3, Distributions.jl version 0.25.108