From 5affef8000bcd71f93259de86257456062b3da29 Mon Sep 17 00:00:00 2001 From: Vincent Gao Date: Thu, 16 Jul 2026 15:58:45 +0200 Subject: [PATCH] fix: correct entropy formulas for Pareto, Geometric, StudentsT, and Dirichlet - Pareto: formula was the exact negation of ln(x_m/a) + 1/a + 1 - Geometric: used log base 2, returning bits instead of nats - StudentsT: scale shifted entropy by -ln(scale) instead of +ln(scale) - Dirichlet: sum of ln_gamma(alpha_i) entered with the wrong sign Existing test assertions for Pareto, Geometric, and Dirichlet encoded the buggy outputs; corrected to closed-form values cross-checked against scipy and extended with more parameter points. Added an entropy test for StudentsT, which had none. --- src/distribution/dirichlet.rs | 28 +++++++++++++++++++--------- src/distribution/geometric.rs | 9 ++++++--- src/distribution/pareto.rs | 17 ++++++++++------- src/distribution/students_t.rs | 19 ++++++++++++++++--- 4 files changed, 51 insertions(+), 22 deletions(-) diff --git a/src/distribution/dirichlet.rs b/src/distribution/dirichlet.rs index 87b6a332..3b7970e0 100644 --- a/src/distribution/dirichlet.rs +++ b/src/distribution/dirichlet.rs @@ -173,11 +173,14 @@ where /// is the `i`th concentration parameter, and `Σ` is the sum from `1` to `K` pub fn entropy(&self) -> Option { let sum = self.alpha_sum(); - let num = self.alpha.iter().fold(0.0, |acc, &x| { - acc + gamma::ln_gamma(x) + (x - 1.0) * gamma::digamma(x) - }); - let entr = - -gamma::ln_gamma(sum) + (sum - self.alpha.len() as f64) * gamma::digamma(sum) - num; + let ln_b = + self.alpha.iter().map(|&x| gamma::ln_gamma(x)).sum::() - gamma::ln_gamma(sum); + let entr = ln_b + (sum - self.alpha.len() as f64) * gamma::digamma(sum) + - self + .alpha + .iter() + .map(|&x| (x - 1.0) * gamma::digamma(x)) + .sum::(); Some(entr) } } @@ -500,16 +503,23 @@ mod tests { #[test] fn test_entropy() { let entropy = |x: Dirichlet<_>| x.entropy().unwrap(); + // cross-checked against scipy.stats.dirichlet(alpha).entropy() test_almost( vector![0.1, 0.3, 0.5, 0.8], - -17.46469081094079, - 1e-30, + -9.318820275187153, + 1e-12, entropy, ); test_almost( vector![0.1, 0.2, 0.3, 0.4], - -21.53881433791513, - 1e-30, + -10.200309764545434, + 1e-12, + entropy, + ); + test_almost( + vector![1.5, 2.0, 3.5, 4.0], + -2.7374675446648483, + 1e-12, entropy, ); } diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs index f471c131..09b431b1 100644 --- a/src/distribution/geometric.rs +++ b/src/distribution/geometric.rs @@ -268,11 +268,11 @@ impl Distribution for Geometric { /// # Formula /// /// ```text - /// (-(1 - p) * log_2(1 - p) - p * log_2(p)) / p + /// (-(1 - p) * ln(1 - p) - p * ln(p)) / p /// ``` fn entropy(&self) -> Option { let inv = 1.0 / self.p; - Some(-inv * (1. - self.p).log(2.0) + (inv - 1.).log(2.0)) + Some(-inv * (1. - self.p).ln() + (inv - 1.).ln()) } /// Returns the skewness of the geometric distribution @@ -396,7 +396,10 @@ mod tests { #[test] fn test_entropy() { let entropy = |x: Geometric| x.entropy().unwrap(); - test_absolute(0.3, 2.937636330768973333333, 1e-14, entropy); + // cross-checked against scipy.stats.geom(p).entropy() (nats) + test_absolute(0.3, 2.0362143401829782, 1e-14, entropy); + test_absolute(0.5, 2.0 * 2f64.ln(), 1e-14, entropy); + test_absolute(0.9, 0.3612033037682758, 1e-14, entropy); test_is_nan(1.0, entropy); } diff --git a/src/distribution/pareto.rs b/src/distribution/pareto.rs index 5b60074b..c9e4d364 100644 --- a/src/distribution/pareto.rs +++ b/src/distribution/pareto.rs @@ -268,12 +268,12 @@ impl Distribution for Pareto { /// # Formula /// /// ```text - /// ln(α/x_m) - 1/α - 1 + /// ln(x_m/α) + 1/α + 1 /// ``` /// /// where `x_m` is the scale and `α` is the shape fn entropy(&self) -> Option { - Some(self.shape.ln() - self.scale.ln() - (1.0 / self.shape) - 1.0) + Some(self.scale.ln() - self.shape.ln() + (1.0 / self.shape) + 1.0) } /// Returns the skewness of the Pareto distribution @@ -423,11 +423,14 @@ mod tests { #[test] fn test_entropy() { let entropy = |x: Pareto| x.entropy().unwrap(); - test_exact(0.1, 0.1, -11.0, entropy); - test_exact(1.0, 1.0, -2.0, entropy); - test_exact(10.0, 10.0, -1.1, entropy); - test_exact(3.0, 1.0, -2.0 - 3f64.ln(), entropy); - test_exact(1.0, 3.0, -4.0/3.0 + 3f64.ln(), entropy); + test_exact(0.1, 0.1, 11.0, entropy); + test_exact(1.0, 1.0, 2.0, entropy); + test_exact(10.0, 10.0, 1.1, entropy); + test_exact(3.0, 1.0, 2.0 + 3f64.ln(), entropy); + test_exact(1.0, 3.0, 4.0 / 3.0 - 3f64.ln(), entropy); + // cross-checked against scipy.stats.pareto(b=shape, scale=scale).entropy() + test_absolute(2.0, 4.5, 0.4112920060058935, 1e-14, entropy); + test_absolute(1.0, 0.5, 3.6931471805599453, 1e-14, entropy); } #[test] diff --git a/src/distribution/students_t.rs b/src/distribution/students_t.rs index 1d0da476..e263c33e 100644 --- a/src/distribution/students_t.rs +++ b/src/distribution/students_t.rs @@ -320,7 +320,7 @@ impl Distribution for StudentsT { /// # Formula /// /// ```text - /// - ln(σ) + (v + 1) / 2 * (ψ((v + 1) / 2) - ψ(v / 2)) + ln(sqrt(v) * B(v / 2, 1 / + /// ln(σ) + (v + 1) / 2 * (ψ((v + 1) / 2) - ψ(v / 2)) + ln(sqrt(v) * B(v / 2, 1 / /// 2)) /// ``` /// @@ -329,8 +329,8 @@ impl Distribution for StudentsT { fn entropy(&self) -> Option { // generalised Student's T is related to normal Student's T by `Y = μ + σ X` // where `X` is distributed as Student's T, plugging into the definition - // of entropy shows scaling affects the entropy by an additive constant `- ln σ` - let shift = -self.scale.ln(); + // of entropy shows scaling affects the entropy by an additive constant `+ ln σ` + let shift = self.scale.ln(); let result = (self.freedom + 1.0) / 2.0 * (gamma::digamma((self.freedom + 1.0) / 2.0) - gamma::digamma(self.freedom / 2.0)) + (self.freedom.sqrt() * beta::beta(self.freedom / 2.0, 0.5)).ln(); @@ -522,6 +522,19 @@ mod tests { test_none(1.0, 1.0, 0.5, |dist| dist.variance()); } + #[test] + fn test_entropy() { + let entropy = |x: StudentsT| x.entropy().unwrap(); + // cross-checked against scipy.stats.t(df, loc, scale).entropy() + test_relative(0.0, 1.0, 3.0, 1.7734775718632907, entropy); + // location does not affect the entropy + test_relative(2.0, 1.0, 3.0, 1.7734775718632907, entropy); + // scale shifts the entropy by +ln(scale) + test_relative(1.0, 2.0, 2.5, 2.5409072511358666, entropy); + test_relative(0.0, 3.0, 5.0, 2.726114961082506, entropy); + test_relative(0.0, 0.5, 10.0, 0.828115312415736, entropy); + } + // TODO: valid skewness tests #[test] fn test_skewness_freedom_lte_3() {