1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.statistics.distribution;
18
19 import java.util.Arrays;
20 import org.apache.commons.math3.stat.StatUtils;
21 import org.apache.commons.math3.stat.inference.GTest;
22 import org.apache.commons.rng.UniformRandomProvider;
23 import org.apache.commons.rng.simple.RandomSource;
24 import org.junit.jupiter.api.Assertions;
25 import org.junit.jupiter.api.Test;
26 import org.junit.jupiter.params.ParameterizedTest;
27 import org.junit.jupiter.params.provider.CsvSource;
28
29
30
31
32
33
34
35
36 class BetaDistributionTest extends BaseContinuousDistributionTest {
37
38 static final double[] ALPHA_BETAS = {0.1, 1, 10, 100, 1000};
39
40 static final double EPSILON = StatUtils.min(ALPHA_BETAS);
41
42 @Override
43 ContinuousDistribution makeDistribution(Object... parameters) {
44 final double alpha = (Double) parameters[0];
45 final double beta = (Double) parameters[1];
46 return BetaDistribution.of(alpha, beta);
47 }
48
49 @Override
50 Object[][] makeInvalidParameters() {
51 return new Object[][] {
52 {0.0, 1.0},
53 {-0.1, 1.0},
54 {0.5, 0.0},
55 {0.5, -0.1}
56 };
57 }
58
59 @Override
60 String[] getParameterNames() {
61 return new String[] {"Alpha", "Beta"};
62 }
63
64 @Override
65 protected double getRelativeTolerance() {
66 return 8e-15;
67 }
68
69
70
71
72
73
74
75 @ParameterizedTest
76 @CsvSource({
77
78 "5.0, 5.0, 0.0001, 1.2595800539968654e-18",
79 "4.0, 5.0, 0.00001, 6.999776002800025e-19",
80 "5.0, 4.0, 0.0001, 5.598600119996539e-19",
81 "6.0, 2.0, 0.001, 6.994000000000028e-18",
82 "2.0, 6.0, 1e-9, 2.0999999930000014e-17",
83 })
84 void testCumulativePrecision(double alpha, double beta, double value, double expected) {
85 final double tolerance = 1e-22;
86 final BetaDistribution dist = BetaDistribution.of(alpha, beta);
87 Assertions.assertEquals(
88 expected,
89 dist.cumulativeProbability(value),
90 tolerance,
91 () -> "cumulative probability not precise at " + value + " for a=" + alpha + " & b=" + beta);
92 }
93
94
95
96
97
98 @ParameterizedTest
99 @CsvSource({
100
101 "5.0, 5.0, 0.9999, 1.2595800539961496e-18",
102 "4.0, 5.0, 0.9999, 5.598600119993397e-19",
103 "5.0, 4.0, 0.99998, 1.1199283217964632e-17",
104 "6.0, 2.0, 0.999999999, 2.0999998742158932e-17",
105 "2.0, 6.0, 0.999, 6.994000000000077e-18",
106 })
107 void testSurvivalPrecision(double alpha, double beta, double value, double expected) {
108 final double tolerance = 1e-22;
109 final BetaDistribution dist = BetaDistribution.of(alpha, beta);
110 Assertions.assertEquals(
111 expected,
112 dist.survivalProbability(value),
113 tolerance,
114 () -> "survival function not precise at " + value + " for a=" + alpha + " & b=" + beta);
115 }
116
117 @ParameterizedTest
118 @CsvSource({
119 "0.5, 3, 0, Infinity",
120 "2, 0.5, 1, Infinity",
121 })
122 void testLogDensityPrecondition(double a, double b, double x, double expected) {
123 final BetaDistribution dist = BetaDistribution.of(a, b);
124 Assertions.assertEquals(expected, dist.density(x));
125 Assertions.assertEquals(Math.log(expected), dist.logDensity(x));
126 }
127
128 @Test
129 void testMomentsSampling() {
130 final UniformRandomProvider rng = RandomSource.XO_SHI_RO_256_PP.create(123456789L);
131 final int numSamples = 1000;
132 for (final double alpha : ALPHA_BETAS) {
133 for (final double beta : ALPHA_BETAS) {
134 final BetaDistribution betaDistribution = BetaDistribution.of(alpha, beta);
135
136 final ContinuousDistribution.Sampler sampler = betaDistribution.createSampler(rng);
137 final double[] observed = TestUtils.sample(numSamples, sampler);
138 Arrays.sort(observed);
139
140 Assertions.assertEquals(betaDistribution.getMean(), StatUtils.mean(observed), EPSILON,
141 () -> String.format("E[Beta(%.2f, %.2f)]", alpha, beta));
142 Assertions.assertEquals(betaDistribution.getVariance(), StatUtils.variance(observed), EPSILON,
143 () -> String.format("Var[Beta(%.2f, %.2f)]", alpha, beta));
144 }
145 }
146 }
147
148 @Test
149 void testGoodnessOfFit() {
150 final UniformRandomProvider rng = RandomSource.XO_SHI_RO_256_PP.create(123456789L);
151
152 final int numSamples = 1000;
153 final double level = 0.01;
154 for (final double alpha : ALPHA_BETAS) {
155 for (final double beta : ALPHA_BETAS) {
156 final BetaDistribution betaDistribution = BetaDistribution.of(alpha, beta);
157
158 final ContinuousDistribution.Sampler sampler = betaDistribution.createSampler(rng);
159 final double[] observed = TestUtils.sample(numSamples, sampler);
160
161 final double gT = gTest(betaDistribution, observed);
162 Assertions.assertFalse(gT < level,
163 () -> String.format("Beta(%s, %s): G goodness-of-fit (%s) test rejected null at alpha = %s",
164 alpha, beta, gT, level));
165 }
166 }
167 }
168
169 private static double gTest(final ContinuousDistribution expectedDistribution, final double[] values) {
170 final int numBins = values.length / 30;
171 final double[] breaks = new double[numBins];
172 for (int b = 0; b < numBins; b++) {
173 breaks[b] = expectedDistribution.inverseCumulativeProbability((double) (b + 1) / numBins);
174 }
175
176 final long[] observed = new long[numBins];
177 for (final double value : values) {
178 int b = Arrays.binarySearch(breaks, value);
179 if (b < 0) {
180 b = -(b + 1);
181 }
182 observed[b]++;
183 }
184
185 final double[] expected = new double[numBins];
186
187
188
189
190 double x0 = 0;
191 for (int b = 0; b < numBins; b++) {
192 final double x1 = breaks[b];
193 expected[b] = expectedDistribution.probability(x0, x1);
194 x0 = x1;
195 }
196
197 return new GTest().gTest(expected, observed);
198 }
199 }