From cd367dc3fb3c8a73361e593a40768f30ffeaa658 Mon Sep 17 00:00:00 2001 From: Mikhail Pak Date: Tue, 17 Oct 2017 21:44:59 +0200 Subject: [PATCH 1/6] fix: formula for the mode of the chi squared distribution Mode cannot be negative --- src/distribution/chi_squared.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/distribution/chi_squared.rs b/src/distribution/chi_squared.rs index d4aeface..a8e41152 100644 --- a/src/distribution/chi_squared.rs +++ b/src/distribution/chi_squared.rs @@ -289,12 +289,12 @@ impl Mode for ChiSquared { /// # Formula /// /// ```ignore - /// k - 2 + /// max(k - 2, 0) /// ``` /// /// where `k` is the degrees of freedom fn mode(&self) -> f64 { - self.g.mode() + self.g.mode().max(0.0) } } @@ -330,8 +330,8 @@ impl Continuous for ChiSquared { #[cfg(test)] mod test { use std::f64; - use statistics::Median; use distribution::ChiSquared; + use statistics::{Median, Mode}; use distribution::internal::*; fn try_create(freedom: f64) -> ChiSquared { From 82685a26910f22086a81be5f667484e3c52ba3e7 Mon Sep 17 00:00:00 2001 From: Mikhail Pak Date: Wed, 18 Oct 2017 09:13:19 +0200 Subject: [PATCH 2/6] fix: chi squared distribution pdf and log pdf --- src/distribution/chi_squared.rs | 84 +++++++++++++++++++++++++++++++-- 1 file changed, 79 insertions(+), 5 deletions(-) diff --git a/src/distribution/chi_squared.rs b/src/distribution/chi_squared.rs index a8e41152..fa5d6efc 100644 --- a/src/distribution/chi_squared.rs +++ b/src/distribution/chi_squared.rs @@ -310,7 +310,11 @@ impl Continuous for ChiSquared { /// /// where `k` is the degrees of freedom and `Γ` is the gamma function fn pdf(&self, x: f64) -> f64 { - self.g.pdf(x) + if x > 0.0 { + self.g.pdf(x) + } else { + 0.0 + } } /// Calculates the log probability density function for the chi-squared @@ -322,7 +326,11 @@ impl Continuous for ChiSquared { /// ln(1 / (2^(k / 2) * Γ(k / 2)) * x^((k / 2) - 1) * e^(-x / 2)) /// ``` fn ln_pdf(&self, x: f64) -> f64 { - self.g.ln_pdf(x) + if x > 0.0 { + self.g.ln_pdf(x) + } else { + -f64::INFINITY + } } } @@ -330,8 +338,8 @@ impl Continuous for ChiSquared { #[cfg(test)] mod test { use std::f64; - use distribution::ChiSquared; use statistics::{Median, Mode}; + use distribution::{ChiSquared, Continuous}; use distribution::internal::*; fn try_create(freedom: f64) -> ChiSquared { @@ -365,10 +373,76 @@ mod test { test_case(3.0, 3.0 - 2.0 / 3.0, |x| x.median()); } + #[test] + fn test_mode() { + test_case(1.0, 0.0, |x| x.mode()); + test_case(2.0, 0.0, |x| x.mode()); + test_case(3.0, 1.0, |x| x.mode()); + test_case(4.5, 2.5, |x| x.mode()); + test_case(10.0, 8.0, |x| x.mode()); + } + + #[test] + fn test_pdf() { + test_case(1.0, 0.0, |x| x.pdf(0.0)); + test_almost(1.0, 1.2000389484301359798, 1e-15, |x| x.pdf(0.1)); + test_almost(1.0, 0.24197072451914334980, 1e-15, |x| x.pdf(1.0)); + test_almost(1.0, 0.010874740337283141714, 1e-15, |x| x.pdf(5.5)); + test_almost(1.0, 4.7000792147504127122e-26, 1e-38, |x| x.pdf(110.1)); + test_case(1.0, 0.0, |x| x.pdf(f64::INFINITY)); + test_case(2.0, 0.0, |x| x.pdf(0.0)); + test_almost(2.0, 0.47561471225035700455, 1e-15, |x| x.pdf(0.1)); + test_almost(2.0, 0.30326532985631671180, 1e-15, |x| x.pdf(1.0)); + test_almost(2.0, 0.031963930603353786351, 1e-15, |x| x.pdf(5.5)); + test_almost(2.0, 6.1810004550085248492e-25, 1e-37, |x| x.pdf(110.1)); + test_case(2.0, 0.0, |x| x.pdf(f64::INFINITY)); + test_case(2.5, 0.0, |x| x.pdf(0.0)); + test_almost(2.5, 0.24812852712543073541, 1e-15, |x| x.pdf(0.1)); + test_almost(2.5, 0.28134822576318228131, 1e-15, |x| x.pdf(1.0)); + test_almost(2.5, 0.045412171451573920401, 1e-15, |x| x.pdf(5.5)); + test_almost(2.5, 1.8574923023527248767e-24, 1e-36, |x| x.pdf(110.1)); + test_case(2.5, 0.0, |x| x.pdf(f64::INFINITY)); + // test_case(f64::INFINITY, 0.0, |x| x.pdf(0.0)); + // test_case(f64::INFINITY, 0.0, |x| x.pdf(0.1)); + // test_case(f64::INFINITY, 0.0, |x| x.pdf(1.0)); + // test_case(f64::INFINITY, 0.0, |x| x.pdf(5.5)); + // test_case(f64::INFINITY, 0.0, |x| x.pdf(110.1)); + // test_case(f64::INFINITY, 0.0, |x| x.pdf(f64::INFINITY)); + } + + #[test] + fn test_ln_pdf() { + test_case(1.0, -f64::INFINITY, |x| x.ln_pdf(0.0)); + test_almost(1.0, 0.18235401329235010023, 1e-15, |x| x.ln_pdf(0.1)); + test_almost(1.0, -1.4189385332046727418, 1e-15, |x| x.ln_pdf(1.0)); + test_almost(1.0, -4.5213125793238853591, 1e-15, |x| x.ln_pdf(5.5)); + test_almost(1.0, -58.319633055068989881, 1e-13, |x| x.ln_pdf(110.1)); + test_case(1.0, -f64::INFINITY, |x| x.ln_pdf(f64::INFINITY)); + test_case(2.0, -f64::INFINITY, |x| x.ln_pdf(0.0)); + test_almost(2.0, -0.74314718055994530942, 1e-15, |x| x.ln_pdf(0.1)); + test_almost(2.0, -1.1931471805599453094, 1e-15, |x| x.ln_pdf(1.0)); + test_almost(2.0, -3.4431471805599453094, 1e-15, |x| x.ln_pdf(5.5)); + test_almost(2.0, -55.743147180559945309, 1e-13, |x| x.ln_pdf(110.1)); + test_case(2.0, -f64::INFINITY, |x| x.ln_pdf(f64::INFINITY)); + test_case(2.5, -f64::INFINITY, |x| x.ln_pdf(0.0)); + test_almost(2.5, -1.3938084125266298963, 1e-15, |x| x.ln_pdf(0.1)); + test_almost(2.5, -1.2681621392781184753, 1e-15, |x| x.ln_pdf(1.0)); + test_almost(2.5, -3.0919751162185121666, 1e-15, |x| x.ln_pdf(5.5)); + test_almost(2.5, -54.642814878345959906, 1e-13, |x| x.ln_pdf(110.1)); + test_case(2.5, -f64::INFINITY, |x| x.ln_pdf(f64::INFINITY)); + // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(0.0)); + // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(0.1)); + // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(1.0)); + // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(5.5)); + // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(110.1)); + // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(f64::INFINITY)); + } + #[test] fn test_continuous() { - // TODO: figure out why this test fails: - //test::check_continuous_distribution(&try_create(1.0), 0.0, 10.0); + // Cannot test for `freedom == 1.0`. The pdf is very steep near 0.0, + // leading to problems with numerical integration + test::check_continuous_distribution(&try_create(1.5), 0.0, 10.0); test::check_continuous_distribution(&try_create(2.0), 0.0, 10.0); test::check_continuous_distribution(&try_create(5.0), 0.0, 50.0); } From 34da91ebd83dea6e4b054b126ce2675162b55a90 Mon Sep 17 00:00:00 2001 From: Mikhail Pak Date: Thu, 19 Oct 2017 12:24:47 +0200 Subject: [PATCH 3/6] test: add mean, variance, stddev, skewness, min, max for chi-squared --- src/distribution/chi_squared.rs | 52 ++++++++++++++++++++++++++++----- 1 file changed, 45 insertions(+), 7 deletions(-) diff --git a/src/distribution/chi_squared.rs b/src/distribution/chi_squared.rs index fa5d6efc..7c7c4006 100644 --- a/src/distribution/chi_squared.rs +++ b/src/distribution/chi_squared.rs @@ -338,7 +338,7 @@ impl Continuous for ChiSquared { #[cfg(test)] mod test { use std::f64; - use statistics::{Median, Mode}; + use statistics::*; use distribution::{ChiSquared, Continuous}; use distribution::internal::*; @@ -365,12 +365,35 @@ mod test { } #[test] - fn test_median() { - test_almost(0.5, 0.0857338820301783264746, 1e-16, |x| x.median()); - test_case(1.0, 1.0 - 2.0 / 3.0, |x| x.median()); - test_case(2.0, 2.0 - 2.0 / 3.0, |x| x.median()); - test_case(2.5, 2.5 - 2.0 / 3.0, |x| x.median()); - test_case(3.0, 3.0 - 2.0 / 3.0, |x| x.median()); + fn test_mean() { + test_case(1.0, 1.0, |x| x.mean()); + test_case(2.1, 2.1, |x| x.mean()); + test_case(3.0, 3.0, |x| x.mean()); + test_case(4.5, 4.5, |x| x.mean()); + } + + #[test] + fn test_variance() { + test_case(1.0, 2.0, |x| x.variance()); + test_case(2.1, 4.2, |x| x.variance()); + test_case(3.0, 6.0, |x| x.variance()); + test_case(4.5, 9.0, |x| x.variance()); + } + + #[test] + fn test_std_dev() { + test_case(1.0, 2f64.sqrt(), |x| x.std_dev()); + test_case(2.1, 4.2f64.sqrt(), |x| x.std_dev()); + test_case(3.0, 6f64.sqrt(), |x| x.std_dev()); + test_case(4.5, 3.0, |x| x.std_dev()); + } + + #[test] + fn test_skewness() { + test_almost(1.0, 8f64.sqrt(), 1e-15, |x| x.skewness()); + test_almost(2.1, (8f64/2.1).sqrt(), 1e-15, |x| x.skewness()); + test_almost(3.0, (8f64/3.0).sqrt(), 1e-15, |x| x.skewness()); + test_almost(4.5, (8f64/4.5).sqrt(), 1e-15, |x| x.skewness()); } #[test] @@ -382,6 +405,21 @@ mod test { test_case(10.0, 8.0, |x| x.mode()); } + #[test] + fn test_median() { + test_almost(0.5, 0.0857338820301783264746, 1e-16, |x| x.median()); + test_case(1.0, 1.0 - 2.0 / 3.0, |x| x.median()); + test_case(2.0, 2.0 - 2.0 / 3.0, |x| x.median()); + test_case(2.5, 2.5 - 2.0 / 3.0, |x| x.median()); + test_case(3.0, 3.0 - 2.0 / 3.0, |x| x.median()); + } + + #[test] + fn test_min_max() { + test_case(1.0, 0.0, |x| x.min()); + test_case(1.0, f64::INFINITY, |x| x.max()); + } + #[test] fn test_pdf() { test_case(1.0, 0.0, |x| x.pdf(0.0)); From 6ca79bc3512c50f2effcea2202d0b1ec93ec41a4 Mon Sep 17 00:00:00 2001 From: Mikhail Pak Date: Thu, 19 Oct 2017 13:50:13 +0200 Subject: [PATCH 4/6] test: add test for the entropy of the chi-squared distribution --- src/distribution/chi_squared.rs | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/distribution/chi_squared.rs b/src/distribution/chi_squared.rs index 7c7c4006..a907b970 100644 --- a/src/distribution/chi_squared.rs +++ b/src/distribution/chi_squared.rs @@ -420,6 +420,16 @@ mod test { test_case(1.0, f64::INFINITY, |x| x.max()); } + #[test] + fn test_entropy() { + test_almost(0.1, -15.760926360123200, 1e-13, |x| x.entropy()); + test_almost(1.0, 0.783757110473934, 1e-15, |x| x.entropy()); + test_almost(2.0, 1.693147180559945, 1e-15, |x| x.entropy()); + test_almost(4.0, 2.270362845461478, 1e-13, |x| x.entropy()); + test_almost(16.0, 3.108818195936093, 1e-13, |x| x.entropy()); + test_almost(100.0, 4.061397128938097, 1e-13, |x| x.entropy()); + } + #[test] fn test_pdf() { test_case(1.0, 0.0, |x| x.pdf(0.0)); From 0f9859e455bd9b61a51dfe9f756b86e8c48eb279 Mon Sep 17 00:00:00 2001 From: Mikhail Pak Date: Thu, 19 Oct 2017 14:48:05 +0200 Subject: [PATCH 5/6] test: add cdf unit test for the chi-squared --- src/distribution/chi_squared.rs | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/src/distribution/chi_squared.rs b/src/distribution/chi_squared.rs index a907b970..b4ecae08 100644 --- a/src/distribution/chi_squared.rs +++ b/src/distribution/chi_squared.rs @@ -339,7 +339,7 @@ impl Continuous for ChiSquared { mod test { use std::f64; use statistics::*; - use distribution::{ChiSquared, Continuous}; + use distribution::{ChiSquared, Univariate, Continuous}; use distribution::internal::*; fn try_create(freedom: f64) -> ChiSquared { @@ -430,6 +430,36 @@ mod test { test_almost(100.0, 4.061397128938097, 1e-13, |x| x.entropy()); } + #[test] + fn test_cdf() { + test_case(0.5, 0.0, |x| x.cdf(0.0)); + test_almost(0.5, 5.165553208304657e-01, 1e-15, |x| x.cdf(0.1)); + test_almost(0.5, 8.464864041916774e-01, 1e-15, |x| x.cdf(1.0)); + test_almost(0.5, 9.995079748755537e-01, 1e-15, |x| x.cdf(10.0)); + test_case(0.5, 1.0, |x| x.cdf(100.0)); + + test_case(1.0, 0.0, |x| x.cdf(0.0)); + test_almost(1.0, 2.481703659541508e-01, 1e-15, |x| x.cdf(0.1)); + test_almost(1.0, 6.826894921370859e-01, 1e-15, |x| x.cdf(1.0)); + test_almost(1.0, 9.984345977419975e-01, 1e-15, |x| x.cdf(10.0)); + test_case(1.0, 1.0, |x| x.cdf(100.0)); + + test_case(3.5, 0.0, |x| x.cdf(0.0)); + test_almost(3.5, 3.184413836670893e-03, 1e-16, |x| x.cdf(0.1)); + test_almost(3.5, 1.355315036184170e-01, 1e-15, |x| x.cdf(1.0)); + test_almost(3.5, 9.719601382856896e-01, 1e-15, |x| x.cdf(10.0)); + test_case(3.5, 1.0, |x| x.cdf(100.0)); + + test_case(100.0, 0.0, |x| x.cdf(0.0)); + test_almost(100.0, 2.181059214078525e-32, 1e-42, |x| x.cdf(10.0)); + test_almost(100.0, 6.953305247616201e-06, 1e-18, |x| x.cdf(50.0)); + test_almost(100.0, 5.188083154720433e-01, 1e-13, |x| x.cdf(100.0)); + + test_almost(400.0, 2.024759014847327e-57, 1e-72, |x| x.cdf(100.0)); + test_almost(400.0, 9.995177872404065e-01, 1e-15, |x| x.cdf(500.0)); + test_case(400.0, 1.0, |x| x.cdf(1000.0)); + } + #[test] fn test_pdf() { test_case(1.0, 0.0, |x| x.pdf(0.0)); From 2df72f2b535ed02e60ee07d3dbcf513fe670c039 Mon Sep 17 00:00:00 2001 From: Mikhail Pak Date: Fri, 20 Oct 2017 11:08:58 +0200 Subject: [PATCH 6/6] fix: use NEG_INFINITY --- src/distribution/chi_squared.rs | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/distribution/chi_squared.rs b/src/distribution/chi_squared.rs index b4ecae08..1d956a1c 100644 --- a/src/distribution/chi_squared.rs +++ b/src/distribution/chi_squared.rs @@ -329,7 +329,7 @@ impl Continuous for ChiSquared { if x > 0.0 { self.g.ln_pdf(x) } else { - -f64::INFINITY + f64::NEG_INFINITY } } } @@ -490,30 +490,30 @@ mod test { #[test] fn test_ln_pdf() { - test_case(1.0, -f64::INFINITY, |x| x.ln_pdf(0.0)); + test_case(1.0, f64::NEG_INFINITY, |x| x.ln_pdf(0.0)); test_almost(1.0, 0.18235401329235010023, 1e-15, |x| x.ln_pdf(0.1)); test_almost(1.0, -1.4189385332046727418, 1e-15, |x| x.ln_pdf(1.0)); test_almost(1.0, -4.5213125793238853591, 1e-15, |x| x.ln_pdf(5.5)); test_almost(1.0, -58.319633055068989881, 1e-13, |x| x.ln_pdf(110.1)); - test_case(1.0, -f64::INFINITY, |x| x.ln_pdf(f64::INFINITY)); - test_case(2.0, -f64::INFINITY, |x| x.ln_pdf(0.0)); + test_case(1.0, f64::NEG_INFINITY, |x| x.ln_pdf(f64::INFINITY)); + test_case(2.0, f64::NEG_INFINITY, |x| x.ln_pdf(0.0)); test_almost(2.0, -0.74314718055994530942, 1e-15, |x| x.ln_pdf(0.1)); test_almost(2.0, -1.1931471805599453094, 1e-15, |x| x.ln_pdf(1.0)); test_almost(2.0, -3.4431471805599453094, 1e-15, |x| x.ln_pdf(5.5)); test_almost(2.0, -55.743147180559945309, 1e-13, |x| x.ln_pdf(110.1)); - test_case(2.0, -f64::INFINITY, |x| x.ln_pdf(f64::INFINITY)); - test_case(2.5, -f64::INFINITY, |x| x.ln_pdf(0.0)); + test_case(2.0, f64::NEG_INFINITY, |x| x.ln_pdf(f64::INFINITY)); + test_case(2.5, f64::NEG_INFINITY, |x| x.ln_pdf(0.0)); test_almost(2.5, -1.3938084125266298963, 1e-15, |x| x.ln_pdf(0.1)); test_almost(2.5, -1.2681621392781184753, 1e-15, |x| x.ln_pdf(1.0)); test_almost(2.5, -3.0919751162185121666, 1e-15, |x| x.ln_pdf(5.5)); test_almost(2.5, -54.642814878345959906, 1e-13, |x| x.ln_pdf(110.1)); - test_case(2.5, -f64::INFINITY, |x| x.ln_pdf(f64::INFINITY)); - // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(0.0)); - // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(0.1)); - // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(1.0)); - // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(5.5)); - // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(110.1)); - // test_case(f64::INFINITY, -f64::INFINITY, |x| x.ln_pdf(f64::INFINITY)); + test_case(2.5, f64::NEG_INFINITY, |x| x.ln_pdf(f64::INFINITY)); + // test_case(f64::INFINITY, f64::NEG_INFINITY, |x| x.ln_pdf(0.0)); + // test_case(f64::INFINITY, f64::NEG_INFINITY, |x| x.ln_pdf(0.1)); + // test_case(f64::INFINITY, f64::NEG_INFINITY, |x| x.ln_pdf(1.0)); + // test_case(f64::INFINITY, f64::NEG_INFINITY, |x| x.ln_pdf(5.5)); + // test_case(f64::INFINITY, f64::NEG_INFINITY, |x| x.ln_pdf(110.1)); + // test_case(f64::INFINITY, f64::NEG_INFINITY, |x| x.ln_pdf(f64::INFINITY)); } #[test]