Fix PERT distribution for when mode is very close to (max - min) / 2
#1311
Conversation
LukasKalbertodt
force-pushed
the
fix-pert-distribution
branch
from
f6161bf
to
b76ed8a
Compare
rand_distr/src/pert.rs
Outdated
| @@ -99,7 +99,7 @@ where | |||
|
|
|||
| let range = max - min; | |||
| let mu = (min + max + shape * mode) / (shape + F::from(2.).unwrap()); | |||
| let v = if mu == mode { | |||
| let v = if (mu - mode).abs() < F::epsilon() { | |||
There was a problem hiding this comment.
Epsilon is not the appropriate value to use here: e.g. if min and max are very small this test will be too permissive. I'd suggest something like 2 * ε / mu2εμ.
There was a problem hiding this comment.
Good point. I've read up on floating point comparisons (in particular this) and force pushed some changes. I'm now using 2 * F::EPSILON * max(mu.abs(), mode.abs()). This is something commonly suggested, i.e. scaling the machine epsilon by the maximum of the two inputs.
There was a problem hiding this comment.
This might still go wrong if one of the values is 0. The linked article talks a lot about how that's quite difficult to deal with. And actually, the popular approx crate treats numbers as equal if they differ less than F::epsilon, at least by default. (To the best of my understanding.) This is to deal with numbers close to 0. So even for ULP or relative epsilon comparison, this absolute epsilon comparison is done first. So actually, thinking about it, my initial commit might be fine for small numbers, but it's main problem was probably large numbers (where F::epsilon is smaller than the gap between two neighboring floats).
Uff, this is tricky :/
There was a problem hiding this comment.
I pushed a new version again. This time, the epsilon is scaled by the largest input. I think with this, everything should be fine. In particular, this deals with "catastrophic cancellation": where the inputs (min, max, mode, shape) are quite large, but due to subtracting them from one another, we keep the quite large rounding error even for mu close to 0. This is also recommended in the linked article:
You’ll need to use an absolute epsilon, whose value might be some small multiple of FLT_EPSILON and the inputs to your calculation.
LukasKalbertodt
force-pushed
the
fix-pert-distribution
branch
3 times, most recently
from
880e1ac
to
ebd1d98
Compare
rand_distr/src/pert.rs
Outdated
| // to the whole calculation. | ||
| // | ||
| // https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/ | ||
| let epsilon = shape.max(min.abs()).max(max.abs()) * F::epsilon() * (F::one() + F::one()); |
There was a problem hiding this comment.
Must you call this complicated expression epsilon? And why does it involve shape anyway?
There is already a special condition for this case, as the general
formula breaks for floats. However, the check uses `==` which is fishy
for floats. This commit instead checks if the difference is smaller than
the machine epsilon.
Without this commit, this returns an error (despite being totally valid
parameters for PERT):
Pert::new(0.0, 0.48258883, 0.24129441)
LukasKalbertodt
force-pushed
the
fix-pert-distribution
branch
from
ebd1d98
to
6f8437e
Compare
I think the problem with this is that But... I've looked at this quite a while now and I'm not 100% sure about anything. It doesn't help that I don't know the semantic meaning of |
There is already a special condition for this case, as the general formula breaks for floats. However, the check uses
==which is fishy for floats. This commit instead checks if the difference is smaller than the machine epsilon.Without this commit, this returns an error (despite being totally valid parameters for PERT):