You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2021/12/24 15:10:20 UTC
[commons-statistics] 01/04: Update to use GammaRatio to compute the moments
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-statistics.git
commit 9ff818307ba59dcc86e8846cc4d26844791f1991
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Fri Dec 24 14:23:53 2021 +0000
Update to use GammaRatio to compute the moments
Added additional test cases added where the simple computation of
gamma(mu + 0.5) / gamma(mu) will overflow the gamma function.
---
.../distribution/NakagamiDistribution.java | 7 +++--
.../distribution/NakagamiDistributionTest.java | 34 ++++++++++++++++++++++
.../distribution/test.nakagami.1.properties | 4 +--
.../distribution/test.nakagami.2.properties | 4 +--
.../distribution/test.nakagami.3.properties | 4 +--
5 files changed, 44 insertions(+), 9 deletions(-)
diff --git a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/NakagamiDistribution.java b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/NakagamiDistribution.java
index 0075495..77f9b66 100644
--- a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/NakagamiDistribution.java
+++ b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/NakagamiDistribution.java
@@ -17,6 +17,7 @@
package org.apache.commons.statistics.distribution;
import org.apache.commons.numbers.gamma.Gamma;
+import org.apache.commons.numbers.gamma.GammaRatio;
import org.apache.commons.numbers.gamma.LogGamma;
import org.apache.commons.numbers.gamma.RegularizedGamma;
@@ -64,9 +65,9 @@ public final class NakagamiDistribution extends AbstractContinuousDistribution {
this.omega = omega;
densityPrefactor = 2.0 * Math.pow(mu, mu) / (Gamma.value(mu) * Math.pow(omega, mu));
logDensityPrefactor = LN_2 + Math.log(mu) * mu - LogGamma.value(mu) - Math.log(omega) * mu;
- mean = Gamma.value(mu + 0.5) / Gamma.value(mu) * Math.sqrt(omega / mu);
- final double v = Gamma.value(mu + 0.5) / Gamma.value(mu);
- variance = omega * (1 - 1 / mu * v * v);
+ final double v = GammaRatio.delta(mu, 0.5);
+ mean = Math.sqrt(omega / mu) / v;
+ variance = omega - mean * mean;
}
/**
diff --git a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/NakagamiDistributionTest.java b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/NakagamiDistributionTest.java
index d21b34a..0438773 100644
--- a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/NakagamiDistributionTest.java
+++ b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/NakagamiDistributionTest.java
@@ -18,6 +18,8 @@ package org.apache.commons.statistics.distribution;
import org.junit.jupiter.api.Assertions;
import org.junit.jupiter.api.Test;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.CsvSource;
/**
* Test cases for {@link NakagamiDistribution}.
@@ -53,6 +55,38 @@ class NakagamiDistributionTest extends BaseContinuousDistributionTest {
//-------------------- Additional test cases -------------------------------
+ /**
+ * Test additional moments.
+ * Includes cases where {@code gamma(mu + 0.5) / gamma(mu)} is not computable
+ * directly due to overflow of the gamma function.
+ */
+ @ParameterizedTest
+ @CsvSource({
+ // Generated using matlab
+ "175, 0.75, 0.86540703592357171, 0.0010706621739778321",
+ "175, 1, 0.99928597029814059, 0.0014275495653037762",
+ "175, 1.25, 1.1172356792742391, 0.0017844369566297202",
+ "175, 3.75, 1.9351089605317091, 0.0053533108698891607",
+ "205.25, 0.75, 0.86549814380218737, 0.00091296307496802065",
+ "205.25, 1, 0.99939117261462862, 0.0012172840999573609",
+ "205.25, 1.25, 1.1173532990397681, 0.0015216051249467011",
+ "205.25, 3.75, 1.9353126839415795, 0.0045648153748401032",
+ "305.25, 0.75, 0.865670838787722, 0.00061399887256183283",
+ "305.25, 1.75, 1.32233404855355, 0.0014326640359776099",
+ "305.25, 3.75, 1.9356988416686078, 0.0030699943628091642",
+ "305.25, 12.75, 3.5692523053388152, 0.010437980833551158",
+ "305.25, 25.25, 5.0228805186490098, 0.020671295376248372",
+ })
+ void testAdditionalMoments(double mu, double omega, double mean, double variance) {
+ // Note:
+ // The relative error of the variance is much greater than the mean.
+ // Accuracy of Matlab data requires verification with another source.
+ // Use a moderate threshold.
+ final DoubleTolerance tolerance = DoubleTolerances.relative(2e-10);
+ final NakagamiDistribution dist = NakagamiDistribution.of(mu, omega);
+ testMoments(dist, mean, variance, tolerance);
+ }
+
@Test
void testExtremeLogDensity() {
// XXX: Verify with more test data from a reference distribution
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.1.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.1.properties
index 7cfb3a6..d86c34c 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.1.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.1.properties
@@ -15,8 +15,8 @@
parameters = 0.5 1.0
# Computed using matlab
-mean = 0.797884560802866
-variance = 0.363380227632418
+mean = 0.79788456080286552
+variance = 0.3633802276324184
lower = 0
upper = Infinity
cdf.points = 0, 1e-3, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.2.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.2.properties
index eacb7ab..bcb8dd9 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.2.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.2.properties
@@ -15,8 +15,8 @@
parameters = 1.5 2.0
# Computed using matlab
-mean = 1.302940031741120
-variance = 0.302347273686450
+mean = 1.3029400317411197
+variance = 0.30234727368644987
lower = 0
upper = Infinity
cdf.points = 0, 1e-3, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.3.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.3.properties
index 21d4db9..3506fbc 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.3.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.nakagami.3.properties
@@ -15,8 +15,8 @@
parameters = 0.3333333333333333333 1.0
# Computed using matlab
-mean = 0.729810132407137
-variance = 0.467377170635876
+mean = 0.72981013240713744
+variance = 0.46737717063587647
lower = 0
upper = Infinity
cdf.points = 0, 1e-3, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2