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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.distribution.continuous;
23
24 import java.io.Serializable;
25
26 import org.hipparchus.analysis.UnivariateFunction;
27 import org.hipparchus.analysis.solvers.UnivariateSolverUtils;
28 import org.hipparchus.distribution.RealDistribution;
29 import org.hipparchus.exception.LocalizedCoreFormats;
30 import org.hipparchus.exception.MathIllegalArgumentException;
31 import org.hipparchus.util.FastMath;
32 import org.hipparchus.util.MathUtils;
33
34 /**
35 * Base class for probability distributions on the reals.
36 * <p>
37 * Default implementations are provided for some of the methods
38 * that do not vary from distribution to distribution.
39 */
40 public abstract class AbstractRealDistribution
41 implements RealDistribution, Serializable {
42
43 /** Default absolute accuracy for inverse cumulative computation. */
44 protected static final double DEFAULT_SOLVER_ABSOLUTE_ACCURACY = 1e-9;
45 /** Serializable version identifier */
46 private static final long serialVersionUID = 20160320L;
47
48 /** Inverse cumulative probability accuracy. */
49 private final double solverAbsoluteAccuracy;
50
51 /** Simple constructor.
52 * @param solverAbsoluteAccuracy the absolute accuracy to use when
53 * computing the inverse cumulative probability.
54 */
55 protected AbstractRealDistribution(double solverAbsoluteAccuracy) {
56 this.solverAbsoluteAccuracy = solverAbsoluteAccuracy;
57 }
58
59 /**
60 * Create a real distribution with default solver absolute accuracy.
61 */
62 protected AbstractRealDistribution() {
63 this.solverAbsoluteAccuracy = DEFAULT_SOLVER_ABSOLUTE_ACCURACY;
64 }
65
66 /**
67 * For a random variable {@code X} whose values are distributed according
68 * to this distribution, this method returns {@code P(x0 < X <= x1)}.
69 *
70 * @param x0 Lower bound (excluded).
71 * @param x1 Upper bound (included).
72 * @return the probability that a random variable with this distribution
73 * takes a value between {@code x0} and {@code x1}, excluding the lower
74 * and including the upper endpoint.
75 * @throws MathIllegalArgumentException if {@code x0 > x1}.
76 * The default implementation uses the identity
77 * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
78 */
79 @Override
80 public double probability(double x0,
81 double x1) throws MathIllegalArgumentException {
82 if (x0 > x1) {
83 throw new MathIllegalArgumentException(LocalizedCoreFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
84 x0, x1, true);
85 }
86 return cumulativeProbability(x1) - cumulativeProbability(x0);
87 }
88
89 /**
90 * {@inheritDoc}
91 *
92 * The default implementation returns
93 * <ul>
94 * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
95 * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
96 * </ul>
97 */
98 @Override
99 public double inverseCumulativeProbability(final double p) throws MathIllegalArgumentException {
100 /*
101 * IMPLEMENTATION NOTES
102 * --------------------
103 * Where applicable, use is made of the one-sided Chebyshev inequality
104 * to bracket the root. This inequality states that
105 * P(X - mu >= k * sig) <= 1 / (1 + k^2),
106 * mu: mean, sig: standard deviation. Equivalently
107 * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
108 * F(mu + k * sig) >= k^2 / (1 + k^2).
109 *
110 * For k = sqrt(p / (1 - p)), we find
111 * F(mu + k * sig) >= p,
112 * and (mu + k * sig) is an upper-bound for the root.
113 *
114 * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
115 * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
116 * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
117 * P(X <= mu - k * sig) <= 1 / (1 + k^2),
118 * F(mu - k * sig) <= 1 / (1 + k^2).
119 *
120 * For k = sqrt((1 - p) / p), we find
121 * F(mu - k * sig) <= p,
122 * and (mu - k * sig) is a lower-bound for the root.
123 *
124 * In cases where the Chebyshev inequality does not apply, geometric
125 * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
126 * the root.
127 */
128
129 MathUtils.checkRangeInclusive(p, 0, 1);
130
131 double lowerBound = getSupportLowerBound();
132 if (p == 0.0) {
133 return lowerBound;
134 }
135
136 double upperBound = getSupportUpperBound();
137 if (p == 1.0) {
138 return upperBound;
139 }
140
141 final double mu = getNumericalMean();
142 final double sig = FastMath.sqrt(getNumericalVariance());
143 final boolean chebyshevApplies;
144 chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
145 Double.isInfinite(sig) || Double.isNaN(sig));
146
147 if (lowerBound == Double.NEGATIVE_INFINITY) {
148 if (chebyshevApplies) {
149 lowerBound = mu - sig * FastMath.sqrt((1. - p) / p);
150 } else {
151 lowerBound = -1.0;
152 while (cumulativeProbability(lowerBound) >= p) {
153 lowerBound *= 2.0;
154 }
155 }
156 }
157
158 if (upperBound == Double.POSITIVE_INFINITY) {
159 if (chebyshevApplies) {
160 upperBound = mu + sig * FastMath.sqrt(p / (1. - p));
161 } else {
162 upperBound = 1.0;
163 while (cumulativeProbability(upperBound) < p) {
164 upperBound *= 2.0;
165 }
166 }
167 }
168
169 final UnivariateFunction toSolve = new UnivariateFunction() {
170 /** {@inheritDoc} */
171 @Override
172 public double value(final double x) {
173 return cumulativeProbability(x) - p;
174 }
175 };
176
177 double x = UnivariateSolverUtils.solve(toSolve,
178 lowerBound,
179 upperBound,
180 getSolverAbsoluteAccuracy());
181
182 if (!isSupportConnected()) {
183 /* Test for plateau. */
184 final double dx = getSolverAbsoluteAccuracy();
185 if (x - dx >= getSupportLowerBound()) {
186 double px = cumulativeProbability(x);
187 if (cumulativeProbability(x - dx) == px) {
188 upperBound = x;
189 while (upperBound - lowerBound > dx) {
190 final double midPoint = 0.5 * (lowerBound + upperBound);
191 if (cumulativeProbability(midPoint) < px) {
192 lowerBound = midPoint;
193 } else {
194 upperBound = midPoint;
195 }
196 }
197 return upperBound;
198 }
199 }
200 }
201 return x;
202 }
203
204 /**
205 * Returns the solver absolute accuracy for inverse cumulative computation.
206 * You can override this method in order to use a Brent solver with an
207 * absolute accuracy different from the default.
208 *
209 * @return the maximum absolute error in inverse cumulative probability estimates
210 */
211 protected double getSolverAbsoluteAccuracy() {
212 return solverAbsoluteAccuracy;
213 }
214
215 /**
216 * {@inheritDoc}
217 * <p>
218 * The default implementation simply computes the logarithm of {@code density(x)}.
219 */
220 @Override
221 public double logDensity(double x) {
222 return FastMath.log(density(x));
223 }
224 }
225