Revise pert

[dhardy]

• Added a CHANGELOG.md entry

#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]

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?

[dhardy]

@newpavlov thoughts on Clippy error (see last comment)?

[MichaelOwenDyer]

Looks good, I think the Clippy warning can be allowed in this case. One tiny suggestion: this is a micro-optimization but we could use shape: Option<F> in PertBuilder so that we have lazy evaluation of the default value via self.shape.unwrap_or_else(|| F::from(4.).unwrap());. But it is not immediately clear to me where that call should happen - the shape is needed in both with_mean and with_mode, and unwrapping the value twice is a bit unclean. So it might be more trouble than its worth. I'll leave it up to you!

[newpavlov]

thoughts on Clippy error (see last comment)?

Personally, I agree with the lint. I expect for T::new to return T or Result<T, Error>. I think we should move the new method to PertBuilder and add Pert::new(min: F, max: F, shape: F, mean: F, mode: F) -> Result<Self, PertError>.

[dhardy]

so that we have lazy evaluation of the default value via self.shape.unwrap_or_else(|| F::from(4.).unwrap());.

That expression should optimise down to a constant.

Pert::new(min: F, max: F, shape: F, mean: F, mode: F) -> Result<Self, PertError>

I sincerely hope not. The builder is supposed to make shape optional and provide an either/or choice between mean and mode.

The alternatives are renaming new → build or having multiple constructors:

• Pert::new_with_mean(min: F, max: F, mean: F)
• Pert::new_with_shape_and_mean(min: F, max: F, shape: F, mean: F)
• Pert::new_with_mode(min: F, max: F, mode: F)
• Pert::new_with_shape_and_mode(min: F, max: F, shape: F, mode: F)

[newpavlov]

The alternatives are renaming new → build or having multiple constructors:

I am not familiar with the most commonly used way to construct Pert, so I don't have a particular opinion here. We could have only Pert::new(min: F, max: F) -> Self and refer to the builder API for everything else.

[dhardy]

We could have only Pert::new(min: F, max: F) -> Self

No; mean or mode is required.

[MichaelOwenDyer]

That expression should optimise down to a constant.

Good point, didn't even consider that.

[newpavlov]

No; mean or mode is required.

Then we need at least two methods new_with_shape and new_with_mean.

[dhardy]

Then we need at least two methods ...

We need precisely the four methods I listed above. We can just use those instead of a builder pattern if preferred.

[dhardy]

Opting to exceptionally allow this builder-constructor (approved by @MichaelOwenDyer). Renaming to build is not neat; using the four methods above is not neat. (And since fn new(..) -> Result<Self, E> is common, it's not like you can really just expect the return value to be Self anyway.)