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 2022/01/03 22:02:45 UTC
[commons-statistics] 01/02: Add second reference data for additional 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 aeeb797aa62e1b7f198070b3f7ea70456dea14f2
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Mon Jan 3 21:57:33 2022 +0000
Add second reference data for additional moments
Update computation of variance to avoid the use of the mean squared.
This avoids loss of precision in (omega/mu) by not having to square the
sqrt(omega/mu).
---
.../distribution/NakagamiDistribution.java | 2 +-
.../distribution/NakagamiDistributionTest.java | 39 +++++++++++++++++++++-
2 files changed, 39 insertions(+), 2 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 6a11c69..87df97a 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
@@ -70,7 +70,7 @@ public final class NakagamiDistribution extends AbstractContinuousDistribution {
logDensityPrefactor = LN_2 + Math.log(mu) * mu - LogGamma.value(mu) - Math.log(omega) * mu;
final double v = GammaRatio.delta(mu, 0.5);
mean = Math.sqrt(omega / mu) / v;
- variance = omega - mean * mean;
+ variance = omega - (omega / mu) / v / v;
}
/**
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 0438773..22c1cec 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
@@ -80,13 +80,50 @@ class NakagamiDistributionTest extends BaseContinuousDistributionTest {
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.
+ // variance = omega - mean^2; omega > 0; x > 0; mean > 0
+ // This computation is subject to cancellation due to subtraction of two large
+ // values to approach a result of zero.
// Use a moderate threshold.
final DoubleTolerance tolerance = DoubleTolerances.relative(2e-10);
final NakagamiDistribution dist = NakagamiDistribution.of(mu, omega);
testMoments(dist, mean, variance, tolerance);
}
+ /**
+ * Repeat test of additional moments with alternative source for the expected result.
+ */
+ @ParameterizedTest
+ @CsvSource({
+ // Generated using 128-bit quad precision implementation using Boost C++:
+ // #include <boost/multiprecision/float128.hpp>
+ // #include <boost/math/special_functions/gamma.hpp>
+ // #define quad boost::multiprecision::float128
+ // T v = boost::math::tgamma_delta_ratio(mu, T(0.5));
+ // T mean = sqrt(omega / mu) / v;
+ // T var = omega - (omega / mu) / v / v;
+ "175, 0.75, 0.865407035923572335404337637742305354, 0.00107066217397678136642741884083229635",
+ "175, 1, 0.999285970298141244170512691211913862, 0.0014275495653023751552365584544430618",
+ "175, 1.25, 1.11723567927423980521693795242933784, 0.00178443695662796894404569806805382725",
+ "175, 3.75, 1.93510896053171023839534780723184735, 0.00535331086988390683213709420416109656",
+ "205.25, 0.75, 0.865498143802251959479795150977083271, 0.000912963074856388060643895128688537674",
+ "205.25, 1, 0.999391172614703197622376095323984551, 0.0012172840998085174141918601715848132",
+ "205.25, 1.25, 1.11735329903985129515900415713529348, 0.00152160512476064676773982521448079983",
+ "205.25, 3.75, 1.93531268394172368161190235734322469, 0.00456481537428194030321947564344316985",
+ "305.25, 0.75, 0.865670838787713729127832304174216151, 0.000613998872576147383115881187594898943",
+ "305.25, 1.75, 1.32233404855353739372707758901129787, 0.00143266403601101056060372277105460371",
+ "305.25, 3.75, 1.93569884166858953645398412102636382, 0.00306999436288073691557940593797382064",
+ "305.25, 12.75, 3.56925230533878138370667203279492999, 0.010437980833794505512969980189112608",
+ "305.25, 25.25, 5.02288051864896241877391197174369638, 0.0206712953767302952315679999823609879",
+ })
+ void testAdditionalMoments2(double mu, double omega, double mean, double variance) {
+ // The mean is within 2 ULP.
+ // The variance is closer than the matlab result but the effect of cancellation
+ // prevents high accuracy.
+ final DoubleTolerance tolerance = DoubleTolerances.relative(1e-12);
+ 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