1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.statistics.distribution;
18
19 import java.util.function.IntUnaryOperator;
20 import org.apache.commons.rng.UniformRandomProvider;
21 import org.apache.commons.rng.sampling.distribution.InverseTransformDiscreteSampler;
22
23 /**
24 * Base class for integer-valued discrete distributions. Default
25 * implementations are provided for some of the methods that do not vary
26 * from distribution to distribution.
27 *
28 * <p>This base class provides a default factory method for creating
29 * a {@linkplain DiscreteDistribution.Sampler sampler instance} that uses the
30 * <a href="https://en.wikipedia.org/wiki/Inverse_transform_sampling">
31 * inversion method</a> for generating random samples that follow the
32 * distribution.
33 *
34 * <p>The class provides functionality to evaluate the probability in a range
35 * using either the cumulative probability or the survival probability.
36 * The survival probability is used if both arguments to
37 * {@link #probability(int, int)} are above the median.
38 * Child classes with a known median can override the default {@link #getMedian()}
39 * method.
40 */
41 abstract class AbstractDiscreteDistribution
42 implements DiscreteDistribution {
43 /** Marker value for no median.
44 * This is a long to be outside the value of any possible int valued median. */
45 private static final long NO_MEDIAN = Long.MIN_VALUE;
46
47 /** Cached value of the median. */
48 private long median = NO_MEDIAN;
49
50 /**
51 * Gets the median. This is used to determine if the arguments to the
52 * {@link #probability(int, int)} function are in the upper or lower domain.
53 *
54 * <p>The default implementation calls {@link #inverseCumulativeProbability(double)}
55 * with a value of 0.5.
56 *
57 * @return the median
58 */
59 int getMedian() {
60 long m = median;
61 if (m == NO_MEDIAN) {
62 m = inverseCumulativeProbability(0.5);
63 median = m;
64 }
65 return (int) m;
66 }
67
68 /** {@inheritDoc} */
69 @Override
70 public double probability(int x0,
71 int x1) {
72 if (x0 > x1) {
73 throw new DistributionException(DistributionException.INVALID_RANGE_LOW_GT_HIGH, x0, x1);
74 }
75 // As per the default interface method handle special cases:
76 // x0 = x1 : return 0
77 // x0 + 1 = x1 : return probability(x1)
78 // Long addition avoids overflow
79 if (x0 + 1L >= x1) {
80 return x0 == x1 ? 0.0 : probability(x1);
81 }
82
83 // Use the survival probability when in the upper domain [3]:
84 //
85 // lower median upper
86 // | | |
87 // 1. |------|
88 // x0 x1
89 // 2. |----------|
90 // x0 x1
91 // 3. |--------|
92 // x0 x1
93
94 final double m = getMedian();
95 if (x0 >= m) {
96 return survivalProbability(x0) - survivalProbability(x1);
97 }
98 return cumulativeProbability(x1) - cumulativeProbability(x0);
99 }
100
101 /**
102 * {@inheritDoc}
103 *
104 * <p>The default implementation returns:
105 * <ul>
106 * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
107 * <li>{@link #getSupportUpperBound()} for {@code p = 1}, or</li>
108 * <li>the result of a binary search between the lower and upper bound using
109 * {@link #cumulativeProbability(int) cumulativeProbability(x)}.
110 * The bounds may be bracketed for efficiency.</li>
111 * </ul>
112 *
113 * @throws IllegalArgumentException if {@code p < 0} or {@code p > 1}
114 */
115 @Override
116 public int inverseCumulativeProbability(double p) {
117 ArgumentUtils.checkProbability(p);
118 return inverseProbability(p, 1 - p, false);
119 }
120
121 /**
122 * {@inheritDoc}
123 *
124 * <p>The default implementation returns:
125 * <ul>
126 * <li>{@link #getSupportLowerBound()} for {@code p = 1},</li>
127 * <li>{@link #getSupportUpperBound()} for {@code p = 0}, or</li>
128 * <li>the result of a binary search between the lower and upper bound using
129 * {@link #survivalProbability(int) survivalProbability(x)}.
130 * The bounds may be bracketed for efficiency.</li>
131 * </ul>
132 *
133 * @throws IllegalArgumentException if {@code p < 0} or {@code p > 1}
134 */
135 @Override
136 public int inverseSurvivalProbability(double p) {
137 ArgumentUtils.checkProbability(p);
138 return inverseProbability(1 - p, p, true);
139 }
140
141 /**
142 * Implementation for the inverse cumulative or survival probability.
143 *
144 * @param p Cumulative probability.
145 * @param q Survival probability.
146 * @param complement Set to true to compute the inverse survival probability
147 * @return the value
148 */
149 private int inverseProbability(double p, double q, boolean complement) {
150
151 int lower = getSupportLowerBound();
152 if (p == 0) {
153 return lower;
154 }
155 int upper = getSupportUpperBound();
156 if (q == 0) {
157 return upper;
158 }
159
160 // The binary search sets the upper value to the mid-point
161 // based on fun(x) >= 0. The upper value is returned.
162 //
163 // Create a function to search for x where the upper bound can be
164 // lowered if:
165 // cdf(x) >= p
166 // sf(x) <= q
167 final IntUnaryOperator fun = complement ?
168 x -> Double.compare(q, survivalProbability(x)) :
169 x -> Double.compare(cumulativeProbability(x), p);
170
171 if (lower == Integer.MIN_VALUE) {
172 if (fun.applyAsInt(lower) >= 0) {
173 return lower;
174 }
175 } else {
176 // this ensures:
177 // cumulativeProbability(lower) < p
178 // survivalProbability(lower) > q
179 // which is important for the solving step
180 lower -= 1;
181 }
182
183 // use the one-sided Chebyshev inequality to narrow the bracket
184 // cf. AbstractContinuousDistribution.inverseCumulativeProbability(double)
185 final double mu = getMean();
186 final double sig = Math.sqrt(getVariance());
187 final boolean chebyshevApplies = Double.isFinite(mu) &&
188 ArgumentUtils.isFiniteStrictlyPositive(sig);
189
190 if (chebyshevApplies) {
191 double tmp = mu - sig * Math.sqrt(q / p);
192 if (tmp > lower) {
193 lower = ((int) Math.ceil(tmp)) - 1;
194 }
195 tmp = mu + sig * Math.sqrt(p / q);
196 if (tmp < upper) {
197 upper = ((int) Math.ceil(tmp)) - 1;
198 }
199 }
200
201 return solveInverseProbability(fun, lower, upper);
202 }
203
204 /**
205 * This is a utility function used by {@link
206 * #inverseProbability(double, double, boolean)}. It assumes
207 * that the inverse probability lies in the bracket {@code
208 * (lower, upper]}. The implementation does simple bisection to find the
209 * smallest {@code x} such that {@code fun(x) >= 0}.
210 *
211 * @param fun Probability function.
212 * @param lowerBound Value satisfying {@code fun(lower) < 0}.
213 * @param upperBound Value satisfying {@code fun(upper) >= 0}.
214 * @return the smallest x
215 */
216 private static int solveInverseProbability(IntUnaryOperator fun,
217 int lowerBound,
218 int upperBound) {
219 // Use long to prevent overflow during computation of the middle
220 long lower = lowerBound;
221 long upper = upperBound;
222 while (lower + 1 < upper) {
223 // Note: Cannot replace division by 2 with a right shift because
224 // (lower + upper) can be negative.
225 final long middle = (lower + upper) / 2;
226 final int pm = fun.applyAsInt((int) middle);
227 if (pm < 0) {
228 lower = middle;
229 } else {
230 upper = middle;
231 }
232 }
233 return (int) upper;
234 }
235
236 /** {@inheritDoc} */
237 @Override
238 public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
239 // Inversion method distribution sampler.
240 return InverseTransformDiscreteSampler.of(rng, this::inverseCumulativeProbability)::sample;
241 }
242 }