rust-random / rand

A Rust library for random number generation.
https://crates.io/crates/rand
Other
1.6k stars 423 forks source link

Revise pert #1452

Open dhardy opened 1 month ago

dhardy commented 1 month ago

1311 show-cases an issue with our PERT implementation. The simplest fix I found was to also use the special-case of v when 2*mode - min - max == 0. A better fix can be found in the R package mc2d, duplicated here.

Calculations are straightforward:

Our original definition of v is:
v = (mu - min) * mode_diff / ((mode - mu) * (max - min))
Substituting in mu = (min + max + shape * mode) / (shape + 2) :
v = (min + max + shape * mode - min*(shape + 2)) / (shape + 2) * (2*mode - min - max) / ((mode*(shape + 2) - min - max - shape*mode) / (shape + 2)) / range
  = (min + max + shape * mode - min*(shape + 2)) / (shape + 2) * (2*mode - min - max) / ((mode*(shape + 2) - min - max - shape*mode) / (shape + 2)) / range
  = (min + max + shape * mode - min*shape - min*2) * (2*mode - min - max) / (2*mode - min - max) / range
  = (max - min + shape * (mode - min)) / range
  = 1 + shape * (mode - min) / range

Now for w:
w = v * (max - mu) / (mu - min);
  = (max - min + shape * (mode - min)) / range * (max - mu) / (mu - min)
  = (max - min + shape * (mode - min)) / range * (max - (min + max + shape * mode) / (shape + 2)) / ((min + max + shape * mode) / (shape + 2) - min)
  = (max - min + shape * (mode - min)) / range * (max * (shape + 2) - min - max - shape * mode) / (shape + 2) / (min + max + shape * mode - min * (shape + 2)) * (shape + 2)
  = (max - min + shape * (mode - min)) / range * (max * shape + max - min - shape * mode) / (max + shape * mode - min * shape - min)
  = (max - min + shape * (mode - min)) * (max - min + shape * (max - mode)) / (max - min + shape * (mode - min)) / range
  = (max - min + shape * (max - mode)) / range
  = 1 + shape * (max - mode) / range

As the R source shows, we can easily calculate mode from mean. Since our argument order (min, max, mode) is already unusual, I replaced the constructor with a builder pattern.

dhardy commented 1 month ago

Clippy says:

error: methods called `new` usually return `Self`
  --> rand_distr/src/pert.rs:91:5
   |
91 | /     pub fn new(min: F, max: F) -> PertBuilder<F> {
92 | |         let shape = F::from(4.0).unwrap();
93 | |         PertBuilder { min, max, shape }
94 | |     }

I don't think Pert::build(min, max).with_mode(mode).unwrap() is better. (And we can't attach a .build() there because it would allow modifying shape after calling .with_mean(mean).) Shall we ignore this?