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 <[email protected]> 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
