1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.legacy.distribution;
18
19 import org.apache.commons.statistics.distribution.DiscreteDistribution;
20 import org.apache.commons.math4.legacy.exception.MathInternalError;
21 import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
22 import org.apache.commons.math4.legacy.exception.OutOfRangeException;
23 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
24 import org.apache.commons.rng.UniformRandomProvider;
25 import org.apache.commons.rng.sampling.distribution.InverseTransformDiscreteSampler;
26 import org.apache.commons.math4.core.jdkmath.JdkMath;
27
28
29
30
31
32
33
34 public abstract class AbstractIntegerDistribution
35 implements DiscreteDistribution {
36
37
38
39
40
41
42
43
44 @Override
45 public double probability(int x0, int x1) throws NumberIsTooLargeException {
46 if (x1 < x0) {
47 throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
48 x0, x1, true);
49 }
50 return cumulativeProbability(x1) - cumulativeProbability(x0);
51 }
52
53
54
55
56
57
58
59
60
61
62
63
64 @Override
65 public int inverseCumulativeProbability(final double p) throws OutOfRangeException {
66 if (p < 0.0 || p > 1.0) {
67 throw new OutOfRangeException(p, 0, 1);
68 }
69
70 int lower = getSupportLowerBound();
71 if (p == 0.0) {
72 return lower;
73 }
74 if (lower == Integer.MIN_VALUE) {
75 if (checkedCumulativeProbability(lower) >= p) {
76 return lower;
77 }
78 } else {
79 lower -= 1;
80
81 }
82
83 int upper = getSupportUpperBound();
84 if (p == 1.0) {
85 return upper;
86 }
87
88
89
90 final double mu = getMean();
91 final double sigma = JdkMath.sqrt(getVariance());
92 final boolean chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
93 Double.isInfinite(sigma) || Double.isNaN(sigma) || sigma == 0.0);
94 if (chebyshevApplies) {
95 double k = JdkMath.sqrt((1.0 - p) / p);
96 double tmp = mu - k * sigma;
97 if (tmp > lower) {
98 lower = ((int) JdkMath.ceil(tmp)) - 1;
99 }
100 k = 1.0 / k;
101 tmp = mu + k * sigma;
102 if (tmp < upper) {
103 upper = ((int) JdkMath.ceil(tmp)) - 1;
104 }
105 }
106
107 return solveInverseCumulativeProbability(p, lower, upper);
108 }
109
110
111
112
113
114
115
116
117
118
119
120
121
122 protected int solveInverseCumulativeProbability(final double p, int lower, int upper) {
123 while (lower + 1 < upper) {
124 int xm = (lower + upper) / 2;
125 if (xm < lower || xm > upper) {
126
127
128
129
130
131 xm = lower + (upper - lower) / 2;
132 }
133
134 double pm = checkedCumulativeProbability(xm);
135 if (pm >= p) {
136 upper = xm;
137 } else {
138 lower = xm;
139 }
140 }
141 return upper;
142 }
143
144
145
146
147
148
149
150
151
152
153
154
155 private double checkedCumulativeProbability(int argument)
156 throws MathInternalError {
157 final double result = cumulativeProbability(argument);
158 if (Double.isNaN(result)) {
159 throw new MathInternalError(LocalizedFormats
160 .DISCRETE_CUMULATIVE_PROBABILITY_RETURNED_NAN, argument);
161 }
162 return result;
163 }
164
165
166
167
168
169
170 @Override
171 public double logProbability(int x) {
172 return JdkMath.log(probability(x));
173 }
174
175
176
177
178
179
180
181
182
183 public static int[] sample(int n,
184 DiscreteDistribution.Sampler sampler) {
185 final int[] samples = new int[n];
186 for (int i = 0; i < n; i++) {
187 samples[i] = sampler.sample();
188 }
189 return samples;
190 }
191
192
193 @Override
194 public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
195
196 return InverseTransformDiscreteSampler.of(rng, this::inverseCumulativeProbability)::sample;
197 }
198 }