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.DoubleBinaryOperator;
20 import java.util.function.DoubleUnaryOperator;
21 import org.apache.commons.numbers.rootfinder.BrentSolver;
22 import org.apache.commons.rng.UniformRandomProvider;
23 import org.apache.commons.rng.sampling.distribution.InverseTransformContinuousSampler;
24
25 /**
26 * Base class for probability distributions on the reals.
27 * Default implementations are provided for some of the methods
28 * that do not vary from distribution to distribution.
29 *
30 * <p>This base class provides a default factory method for creating
31 * a {@linkplain ContinuousDistribution.Sampler sampler instance} that uses the
32 * <a href="https://en.wikipedia.org/wiki/Inverse_transform_sampling">
33 * inversion method</a> for generating random samples that follow the
34 * distribution.
35 *
36 * <p>The class provides functionality to evaluate the probability in a range
37 * using either the cumulative probability or the survival probability.
38 * The survival probability is used if both arguments to
39 * {@link #probability(double, double)} are above the median.
40 * Child classes with a known median can override the default {@link #getMedian()}
41 * method.
42 */
43 abstract class AbstractContinuousDistribution
44 implements ContinuousDistribution {
45
46 // Notes on the inverse probability implementation:
47 //
48 // The Brent solver does not allow a stopping criteria for the proximity
49 // to the root; it uses equality to zero within 1 ULP. The search is
50 // iterated until there is a small difference between the upper
51 // and lower bracket of the root, expressed as a combination of relative
52 // and absolute thresholds.
53
54 /** BrentSolver relative accuracy.
55 * This is used with {@code tol = 2 * relEps * abs(b) + absEps} so the minimum
56 * non-zero value with an effect is half of machine epsilon (2^-53). */
57 private static final double SOLVER_RELATIVE_ACCURACY = 0x1.0p-53;
58 /** BrentSolver absolute accuracy.
59 * This is used with {@code tol = 2 * relEps * abs(b) + absEps} so set to MIN_VALUE
60 * so that when the relative epsilon has no effect (as b is too small) the tolerance
61 * is at least 1 ULP for sub-normal numbers. */
62 private static final double SOLVER_ABSOLUTE_ACCURACY = Double.MIN_VALUE;
63 /** BrentSolver function value accuracy.
64 * Determines if the Brent solver performs a search. It is not used during the search.
65 * Set to a very low value to search using Brent's method unless
66 * the starting point is correct, or within 1 ULP for sub-normal probabilities. */
67 private static final double SOLVER_FUNCTION_VALUE_ACCURACY = Double.MIN_VALUE;
68
69 /** Cached value of the median. */
70 private double median = Double.NaN;
71
72 /**
73 * Gets the median. This is used to determine if the arguments to the
74 * {@link #probability(double, double)} function are in the upper or lower domain.
75 *
76 * <p>The default implementation calls {@link #inverseCumulativeProbability(double)}
77 * with a value of 0.5.
78 *
79 * @return the median
80 */
81 double getMedian() {
82 double m = median;
83 if (Double.isNaN(m)) {
84 m = inverseCumulativeProbability(0.5);
85 median = m;
86 }
87 return m;
88 }
89
90 /** {@inheritDoc} */
91 @Override
92 public double probability(double x0,
93 double x1) {
94 if (x0 > x1) {
95 throw new DistributionException(DistributionException.INVALID_RANGE_LOW_GT_HIGH, x0, x1);
96 }
97 // Use the survival probability when in the upper domain [3]:
98 //
99 // lower median upper
100 // | | |
101 // 1. |------|
102 // x0 x1
103 // 2. |----------|
104 // x0 x1
105 // 3. |--------|
106 // x0 x1
107
108 final double m = getMedian();
109 if (x0 >= m) {
110 return survivalProbability(x0) - survivalProbability(x1);
111 }
112 return cumulativeProbability(x1) - cumulativeProbability(x0);
113 }
114
115 /**
116 * {@inheritDoc}
117 *
118 * <p>The default implementation returns:
119 * <ul>
120 * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
121 * <li>{@link #getSupportUpperBound()} for {@code p = 1}, or</li>
122 * <li>the result of a search for a root between the lower and upper bound using
123 * {@link #cumulativeProbability(double) cumulativeProbability(x) - p}.
124 * The bounds may be bracketed for efficiency.</li>
125 * </ul>
126 *
127 * @throws IllegalArgumentException if {@code p < 0} or {@code p > 1}
128 */
129 @Override
130 public double inverseCumulativeProbability(double p) {
131 ArgumentUtils.checkProbability(p);
132 return inverseProbability(p, 1 - p, false);
133 }
134
135 /**
136 * {@inheritDoc}
137 *
138 * <p>The default implementation returns:
139 * <ul>
140 * <li>{@link #getSupportLowerBound()} for {@code p = 1},</li>
141 * <li>{@link #getSupportUpperBound()} for {@code p = 0}, or</li>
142 * <li>the result of a search for a root between the lower and upper bound using
143 * {@link #survivalProbability(double) survivalProbability(x) - p}.
144 * The bounds may be bracketed for efficiency.</li>
145 * </ul>
146 *
147 * @throws IllegalArgumentException if {@code p < 0} or {@code p > 1}
148 */
149 @Override
150 public double inverseSurvivalProbability(double p) {
151 ArgumentUtils.checkProbability(p);
152 return inverseProbability(1 - p, p, true);
153 }
154
155 /**
156 * Implementation for the inverse cumulative or survival probability.
157 *
158 * @param p Cumulative probability.
159 * @param q Survival probability.
160 * @param complement Set to true to compute the inverse survival probability
161 * @return the value
162 */
163 private double inverseProbability(final double p, final double q, boolean complement) {
164 /* IMPLEMENTATION NOTES
165 * --------------------
166 * Where applicable, use is made of the one-sided Chebyshev inequality
167 * to bracket the root. This inequality states that
168 * P(X - mu >= k * sig) <= 1 / (1 + k^2),
169 * mu: mean, sig: standard deviation. Equivalently
170 * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
171 * F(mu + k * sig) >= k^2 / (1 + k^2).
172 *
173 * For k = sqrt(p / (1 - p)), we find
174 * F(mu + k * sig) >= p,
175 * and (mu + k * sig) is an upper-bound for the root.
176 *
177 * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
178 * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
179 * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
180 * P(X <= mu - k * sig) <= 1 / (1 + k^2),
181 * F(mu - k * sig) <= 1 / (1 + k^2).
182 *
183 * For k = sqrt((1 - p) / p), we find
184 * F(mu - k * sig) <= p,
185 * and (mu - k * sig) is a lower-bound for the root.
186 *
187 * In cases where the Chebyshev inequality does not apply, geometric
188 * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
189 * the root.
190 *
191 * In the case of the survival probability the bracket can be set using the same
192 * bound given that the argument p = 1 - q, with q the survival probability.
193 */
194
195 double lowerBound = getSupportLowerBound();
196 if (p == 0) {
197 return lowerBound;
198 }
199 double upperBound = getSupportUpperBound();
200 if (q == 0) {
201 return upperBound;
202 }
203
204 final double mu = getMean();
205 final double sig = Math.sqrt(getVariance());
206 final boolean chebyshevApplies = Double.isFinite(mu) &&
207 ArgumentUtils.isFiniteStrictlyPositive(sig);
208
209 if (lowerBound == Double.NEGATIVE_INFINITY) {
210 lowerBound = createFiniteLowerBound(p, q, complement, upperBound, mu, sig, chebyshevApplies);
211 }
212
213 if (upperBound == Double.POSITIVE_INFINITY) {
214 upperBound = createFiniteUpperBound(p, q, complement, lowerBound, mu, sig, chebyshevApplies);
215 }
216
217 // Here the bracket [lower, upper] uses finite values. If the support
218 // is infinite the bracket can truncate the distribution and the target
219 // probability can be outside the range of [lower, upper].
220 if (upperBound == Double.MAX_VALUE) {
221 if (complement) {
222 if (survivalProbability(upperBound) > q) {
223 return getSupportUpperBound();
224 }
225 } else if (cumulativeProbability(upperBound) < p) {
226 return getSupportUpperBound();
227 }
228 }
229 if (lowerBound == -Double.MAX_VALUE) {
230 if (complement) {
231 if (survivalProbability(lowerBound) < q) {
232 return getSupportLowerBound();
233 }
234 } else if (cumulativeProbability(lowerBound) > p) {
235 return getSupportLowerBound();
236 }
237 }
238
239 final DoubleUnaryOperator fun = complement ?
240 arg -> survivalProbability(arg) - q :
241 arg -> cumulativeProbability(arg) - p;
242 // Note the initial value is robust to overflow.
243 // Do not use 0.5 * (lowerBound + upperBound).
244 final double x = new BrentSolver(SOLVER_RELATIVE_ACCURACY,
245 SOLVER_ABSOLUTE_ACCURACY,
246 SOLVER_FUNCTION_VALUE_ACCURACY)
247 .findRoot(fun,
248 lowerBound,
249 lowerBound + 0.5 * (upperBound - lowerBound),
250 upperBound);
251
252 if (!isSupportConnected()) {
253 return searchPlateau(complement, lowerBound, x);
254 }
255 return x;
256 }
257
258 /**
259 * Create a finite lower bound. Assumes the current lower bound is negative infinity.
260 *
261 * @param p Cumulative probability.
262 * @param q Survival probability.
263 * @param complement Set to true to compute the inverse survival probability
264 * @param upperBound Current upper bound
265 * @param mu Mean
266 * @param sig Standard deviation
267 * @param chebyshevApplies True if the Chebyshev inequality applies (mean is finite and {@code sig > 0}}
268 * @return the finite lower bound
269 */
270 private double createFiniteLowerBound(final double p, final double q, boolean complement,
271 double upperBound, final double mu, final double sig, final boolean chebyshevApplies) {
272 double lowerBound;
273 if (chebyshevApplies) {
274 lowerBound = mu - sig * Math.sqrt(q / p);
275 } else {
276 lowerBound = Double.NEGATIVE_INFINITY;
277 }
278 // Bound may have been set as infinite
279 if (lowerBound == Double.NEGATIVE_INFINITY) {
280 lowerBound = Math.min(-1, upperBound);
281 if (complement) {
282 while (survivalProbability(lowerBound) < q) {
283 lowerBound *= 2;
284 }
285 } else {
286 while (cumulativeProbability(lowerBound) >= p) {
287 lowerBound *= 2;
288 }
289 }
290 // Ensure finite
291 lowerBound = Math.max(lowerBound, -Double.MAX_VALUE);
292 }
293 return lowerBound;
294 }
295
296 /**
297 * Create a finite upper bound. Assumes the current upper bound is positive infinity.
298 *
299 * @param p Cumulative probability.
300 * @param q Survival probability.
301 * @param complement Set to true to compute the inverse survival probability
302 * @param lowerBound Current lower bound
303 * @param mu Mean
304 * @param sig Standard deviation
305 * @param chebyshevApplies True if the Chebyshev inequality applies (mean is finite and {@code sig > 0}}
306 * @return the finite lower bound
307 */
308 private double createFiniteUpperBound(final double p, final double q, boolean complement,
309 double lowerBound, final double mu, final double sig, final boolean chebyshevApplies) {
310 double upperBound;
311 if (chebyshevApplies) {
312 upperBound = mu + sig * Math.sqrt(p / q);
313 } else {
314 upperBound = Double.POSITIVE_INFINITY;
315 }
316 // Bound may have been set as infinite
317 if (upperBound == Double.POSITIVE_INFINITY) {
318 upperBound = Math.max(1, lowerBound);
319 if (complement) {
320 while (survivalProbability(upperBound) >= q) {
321 upperBound *= 2;
322 }
323 } else {
324 while (cumulativeProbability(upperBound) < p) {
325 upperBound *= 2;
326 }
327 }
328 // Ensure finite
329 upperBound = Math.min(upperBound, Double.MAX_VALUE);
330 }
331 return upperBound;
332 }
333
334 /**
335 * Indicates whether the support is connected, i.e. whether all values between the
336 * lower and upper bound of the support are included in the support.
337 *
338 * <p>This method is used in the default implementation of the inverse cumulative and
339 * survival probability functions.
340 *
341 * <p>The default value is true which assumes the cdf and sf have no plateau regions
342 * where the same probability value is returned for a large range of x.
343 * Override this method if there are gaps in the support of the cdf and sf.
344 *
345 * <p>If false then the inverse will perform an additional step to ensure that the
346 * lower-bound of the interval on which the cdf is constant should be returned. This
347 * will search from the initial point x downwards if a smaller value also has the same
348 * cumulative (survival) probability.
349 *
350 * <p>Any plateau with a width in x smaller than the inverse absolute accuracy will
351 * not be searched.
352 *
353 * <p>Note: This method was public in commons math. It has been reduced to package private
354 * in commons statistics as it is an implementation detail.
355 *
356 * @return whether the support is connected.
357 * @see <a href="https://issues.apache.org/jira/browse/MATH-699">MATH-699</a>
358 */
359 boolean isSupportConnected() {
360 return true;
361 }
362
363 /**
364 * Test the probability function for a plateau at the point x. If detected
365 * search the plateau for the lowest point y such that
366 * {@code inf{y in R | P(y) == P(x)}}.
367 *
368 * <p>This function is used when the distribution support is not connected
369 * to satisfy the inverse probability requirements of {@link ContinuousDistribution}
370 * on the returned value.
371 *
372 * @param complement Set to true to search the survival probability.
373 * @param lower Lower bound used to limit the search downwards.
374 * @param x Current value.
375 * @return the infimum y
376 */
377 private double searchPlateau(boolean complement, double lower, final double x) {
378 // Test for plateau. Lower the value x if the probability is the same.
379 // Ensure the step is robust to the solver accuracy being less
380 // than 1 ulp of x (e.g. dx=0 will infinite loop)
381 final double dx = Math.max(SOLVER_ABSOLUTE_ACCURACY, Math.ulp(x));
382 if (x - dx >= lower) {
383 final DoubleUnaryOperator fun = complement ?
384 this::survivalProbability :
385 this::cumulativeProbability;
386 final double px = fun.applyAsDouble(x);
387 if (fun.applyAsDouble(x - dx) == px) {
388 double upperBound = x;
389 double lowerBound = lower;
390 // Bisection search
391 // Require cdf(x) < px and sf(x) > px to move the lower bound
392 // to the midpoint.
393 final DoubleBinaryOperator cmp = complement ?
394 (a, b) -> a > b ? -1 : 1 :
395 (a, b) -> a < b ? -1 : 1;
396 while (upperBound - lowerBound > dx) {
397 final double midPoint = 0.5 * (lowerBound + upperBound);
398 if (cmp.applyAsDouble(fun.applyAsDouble(midPoint), px) < 0) {
399 lowerBound = midPoint;
400 } else {
401 upperBound = midPoint;
402 }
403 }
404 return upperBound;
405 }
406 }
407 return x;
408 }
409
410 /** {@inheritDoc} */
411 @Override
412 public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
413 // Inversion method distribution sampler.
414 return InverseTransformContinuousSampler.of(rng, this::inverseCumulativeProbability)::sample;
415 }
416 }