fix: correct entropy formulas for Pareto, Geometric, StudentsT, and Dirichlet

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<f64> {
         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::<f64>() - 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::<f64>();
         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,
         );
     }

src/distribution/geometric.rs

@@ -268,11 +268,11 @@ impl Distribution<f64> 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<f64> {
         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);
     }
 

src/distribution/pareto.rs

@@ -268,12 +268,12 @@ impl Distribution<f64> 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<f64> {
-        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]

src/distribution/students_t.rs

@@ -320,7 +320,7 @@ impl Distribution<f64> 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<f64> for StudentsT {
     fn entropy(&self) -> Option<f64> {
         // 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() {