New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Removes gamma special cases #104
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Looks good, thanks for doing this! Just one minor nit
This looks good to me, if you're done I can go ahead and merge it in lmk |
…ions fix: homogenised BadParam error checking by replacing match statements with if-elses
I disallowed special parameter cases for other distributions too, the reasons being the same discussed with gamma. |
Looking back at the original issue I filed I should have probably first initiated the discussion with separate issues about those functions. Forgot about that, so I hope you don't mind. |
test_case(1.0, 0.1, 0.036787944117144234201693506390001264039984687455876246, |x| x.pdf(10.0)); | ||
test_case(1.0, 1.0, 0.36787944117144232159552377016146086744581113103176804, |x| x.pdf(1.0)); | ||
test_case(1.0, 1.0, 0.000045399929762484851535591515560550610237918088866564953, |x| x.pdf(10.0)); | ||
test_almost(1.0, 0.1, 0.090483741803595961836995913651194571475319347018875963, 1e-10, |x| x.pdf(1.0)); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm curious as to why these changed and is 1e-10
the strictest bound on them?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I can't explain why those are removed, must be an error on my part.
I'm not sure what a good bound is or how the values in the tests were obtained. Is there a rule to follow on this?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sorry for the confusion, I only saw the four lines outlined above and assumed I'd mistakenly removed tests.
In general testing exact values in floating point is a bad idea. I opened an issue #108 explaining my rationale.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
A lot of these were either ported directly from Math.NET or generated by comparing to other gamma distribution calculators if there were no comparable unit tests in the Math.NET library.
I definitely agree with #108 and as part of that issue it might be worth it to determine what a good bound should be. Anecdotally 1e-14
is probably as much precision as most people will need although there are definitely unit tests that don't meet that criteria (I know there's a couple that are only 1e-12
).
These are probably fine as is for now I was just asking out of curiousity
@@ -308,14 +309,10 @@ impl Continuous<f64, f64> for Weibull { | |||
fn ln_pdf(&self, x: f64) -> f64 { | |||
if x < 0.0 { | |||
f64::NEG_INFINITY | |||
} else if x == 0.0 && self.shape == 1.0 { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
can we add these two paths back in, the 0.0 -
is probably unnecessary but the way they were written is to prevent unnecessary work or divisions
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I've also removed other specialisations, do you want me to put those back in too?
I'll have to profile it. It gets rid of divisions but more logarithms have to be evaluated. An additional division could be left out by writing x.ln() - self.scale.ln()
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't see any other ones that were removed outside the ones that are no longer relevant due to removals of f64::INFINITY
but good point about ln
. Do you mind running a quick benchmark? Should be pretty simple using https://github.com/bluss/bencher/, here's a toy example https://github.com/boxtown/counter/blob/master/benches/bench.rs from another repo
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In the formulas of pdf
and ln_pdf
for gamma I removed the specialisation shape = 1
.
My preference is towards predictability. Trying to manage special cases leads to uglier code with hard to predict control flow. As such tricks accumulate it gets increasingly hard to understand what's going on, with the interactions with caches and branch prediction. On some machines it may become quicker, on some slower.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
running 4 tests
test weibull_new_x_equals_0 ... bench: 4 ns/iter (+/- 0)
test weibull_new_x_equals_2 ... bench: 110 ns/iter (+/- 4)
test weibull_old_x_equals_0 ... bench: 25 ns/iter (+/- 0)
test weibull_old_x_equals_2 ... bench: 117 ns/iter (+/- 8)
Here's the bench if you want to run it.
#[macro_use]
extern crate bencher;
extern crate statrs;
use bencher::Bencher;
use std::f64::*;
use statrs::distribution::{Weibull,Continuous};
struct WeibullNew(Weibull);
impl WeibullNew {
fn ln_pdf(&self, x: f64) -> f64 {
if x < 0.0 {
NEG_INFINITY
} else if x == INFINITY {
NEG_INFINITY
} else {
self.0.pdf(x).ln()
}
}
}
fn weibull_old_x_equals_0(bench: &mut Bencher) {
let w = Weibull::new(1.0, 0.5).unwrap();
bench.iter(|| {
w.ln_pdf(0.0);
})
}
fn weibull_new_x_equals_0(bench: &mut Bencher) {
let w = WeibullNew(Weibull::new(1.0, 0.5).unwrap());
bench.iter(|| {
w.ln_pdf(0.0);
})
}
fn weibull_old_x_equals_2(bench: &mut Bencher) {
let w = Weibull::new(1.0, 0.5).unwrap();
bench.iter(|| {
w.ln_pdf(2.0);
})
}
fn weibull_new_x_equals_2(bench: &mut Bencher) {
let w = WeibullNew(Weibull::new(1.0, 0.5).unwrap());
bench.iter(|| {
w.ln_pdf(2.0);
})
}
benchmark_group!(first, weibull_old_x_equals_0, weibull_new_x_equals_0);
benchmark_group!(second, weibull_old_x_equals_2, weibull_new_x_equals_2);
benchmark_main!(first,second);
Splitting into multiple pr, one for each function. |
Removed also the tests which instantiated gamma with infinity as a parameter.
Tests are still failing, but don't know what the motivation is behind those
numbers.
If it doesn't matter too much I'll just update the testvalues.