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 * http://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.math4.legacy.distribution;
18
19 import org.apache.commons.statistics.distribution.ContinuousDistribution;
20 import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
21 import org.apache.commons.math4.legacy.analysis.solvers.UnivariateSolverUtils;
22 import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
23 import org.apache.commons.math4.legacy.exception.OutOfRangeException;
24 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
25 import org.apache.commons.rng.UniformRandomProvider;
26 import org.apache.commons.rng.sampling.distribution.InverseTransformContinuousSampler;
27 import org.apache.commons.math4.core.jdkmath.JdkMath;
28
29 /**
30 * Base class for probability distributions on the reals.
31 * Default implementations are provided for some of the methods
32 * that do not vary from distribution to distribution.
33 *
34 * <p>
35 * This base class provides a default factory method for creating
36 * a {@link org.apache.commons.statistics.distribution.ContinuousDistribution.Sampler
37 * sampler instance} that uses the
38 * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
39 * inversion method</a> for generating random samples that follow the
40 * distribution.
41 * </p>
42 *
43 * @since 3.0
44 */
45 public abstract class AbstractRealDistribution
46 implements ContinuousDistribution {
47 /** Default absolute accuracy for inverse cumulative computation. */
48 public static final double SOLVER_DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
49
50 /**
51 * For a random variable {@code X} whose values are distributed according
52 * to this distribution, this method returns {@code P(x0 < X <= x1)}.
53 *
54 * @param x0 Lower bound (excluded).
55 * @param x1 Upper bound (included).
56 * @return the probability that a random variable with this distribution
57 * takes a value between {@code x0} and {@code x1}, excluding the lower
58 * and including the upper endpoint.
59 * @throws NumberIsTooLargeException if {@code x0 > x1}.
60 *
61 * The default implementation uses the identity
62 * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
63 *
64 * @since 3.1
65 */
66 @Override
67 public double probability(double x0,
68 double x1) {
69 if (x0 > x1) {
70 throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
71 x0, x1, true);
72 }
73 return cumulativeProbability(x1) - cumulativeProbability(x0);
74 }
75
76 /**
77 * {@inheritDoc}
78 *
79 * The default implementation returns
80 * <ul>
81 * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
82 * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
83 * </ul>
84 */
85 @Override
86 public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
87 /*
88 * IMPLEMENTATION NOTES
89 * --------------------
90 * Where applicable, use is made of the one-sided Chebyshev inequality
91 * to bracket the root. This inequality states that
92 * P(X - mu >= k * sig) <= 1 / (1 + k^2),
93 * mu: mean, sig: standard deviation. Equivalently
94 * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
95 * F(mu + k * sig) >= k^2 / (1 + k^2).
96 *
97 * For k = sqrt(p / (1 - p)), we find
98 * F(mu + k * sig) >= p,
99 * and (mu + k * sig) is an upper-bound for the root.
100 *
101 * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
102 * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
103 * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
104 * P(X <= mu - k * sig) <= 1 / (1 + k^2),
105 * F(mu - k * sig) <= 1 / (1 + k^2).
106 *
107 * For k = sqrt((1 - p) / p), we find
108 * F(mu - k * sig) <= p,
109 * and (mu - k * sig) is a lower-bound for the root.
110 *
111 * In cases where the Chebyshev inequality does not apply, geometric
112 * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
113 * the root.
114 */
115 if (p < 0.0 || p > 1.0) {
116 throw new OutOfRangeException(p, 0, 1);
117 }
118
119 double lowerBound = getSupportLowerBound();
120 if (p == 0.0) {
121 return lowerBound;
122 }
123
124 double upperBound = getSupportUpperBound();
125 if (p == 1.0) {
126 return upperBound;
127 }
128
129 final double mu = getMean();
130 final double sig = JdkMath.sqrt(getVariance());
131 final boolean chebyshevApplies;
132 chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
133 Double.isInfinite(sig) || Double.isNaN(sig));
134
135 if (lowerBound == Double.NEGATIVE_INFINITY) {
136 if (chebyshevApplies) {
137 lowerBound = mu - sig * JdkMath.sqrt((1. - p) / p);
138 } else {
139 lowerBound = -1.0;
140 while (cumulativeProbability(lowerBound) >= p) {
141 lowerBound *= 2.0;
142 }
143 }
144 }
145
146 if (upperBound == Double.POSITIVE_INFINITY) {
147 if (chebyshevApplies) {
148 upperBound = mu + sig * JdkMath.sqrt(p / (1. - p));
149 } else {
150 upperBound = 1.0;
151 while (cumulativeProbability(upperBound) < p) {
152 upperBound *= 2.0;
153 }
154 }
155 }
156
157 final UnivariateFunction toSolve = new UnivariateFunction() {
158 /** {@inheritDoc} */
159 @Override
160 public double value(final double x) {
161 return cumulativeProbability(x) - p;
162 }
163 };
164
165 return UnivariateSolverUtils.solve(toSolve,
166 lowerBound,
167 upperBound,
168 getSolverAbsoluteAccuracy());
169 }
170
171 /**
172 * Returns the solver absolute accuracy for inverse cumulative computation.
173 * You can override this method in order to use a Brent solver with an
174 * absolute accuracy different from the default.
175 *
176 * @return the maximum absolute error in inverse cumulative probability estimates
177 */
178 protected double getSolverAbsoluteAccuracy() {
179 return SOLVER_DEFAULT_ABSOLUTE_ACCURACY;
180 }
181
182 /**
183 * {@inheritDoc}
184 * <p>
185 * The default implementation simply computes the logarithm of {@code density(x)}.
186 */
187 @Override
188 public double logDensity(double x) {
189 return JdkMath.log(density(x));
190 }
191
192 /**
193 * Utility function for allocating an array and filling it with {@code n}
194 * samples generated by the given {@code sampler}.
195 *
196 * @param n Number of samples.
197 * @param sampler Sampler.
198 * @return an array of size {@code n}.
199 */
200 public static double[] sample(int n,
201 ContinuousDistribution.Sampler sampler) {
202 final double[] samples = new double[n];
203 for (int i = 0; i < n; i++) {
204 samples[i] = sampler.sample();
205 }
206 return samples;
207 }
208
209 /**{@inheritDoc} */
210 @Override
211 public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
212 // Inversion method distribution sampler.
213 return InverseTransformContinuousSampler.of(rng, this::inverseCumulativeProbability)::sample;
214 }
215 }