001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.distribution;
018
019 import java.io.Serializable;
020
021 import org.apache.commons.math3.analysis.UnivariateFunction;
022 import org.apache.commons.math3.analysis.solvers.UnivariateSolverUtils;
023 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
024 import org.apache.commons.math3.exception.NumberIsTooLargeException;
025 import org.apache.commons.math3.exception.OutOfRangeException;
026 import org.apache.commons.math3.exception.util.LocalizedFormats;
027 import org.apache.commons.math3.random.RandomGenerator;
028 import org.apache.commons.math3.random.RandomDataImpl;
029 import org.apache.commons.math3.util.FastMath;
030
031 /**
032 * Base class for probability distributions on the reals.
033 * Default implementations are provided for some of the methods
034 * that do not vary from distribution to distribution.
035 *
036 * @version $Id: AbstractRealDistribution.java 1422195 2012-12-15 06:45:18Z psteitz $
037 * @since 3.0
038 */
039 public abstract class AbstractRealDistribution
040 implements RealDistribution, Serializable {
041 /** Default accuracy. */
042 public static final double SOLVER_DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
043 /** Serializable version identifier */
044 private static final long serialVersionUID = -38038050983108802L;
045 /**
046 * RandomData instance used to generate samples from the distribution.
047 * @deprecated As of 3.1, to be removed in 4.0. Please use the
048 * {@link #random} instance variable instead.
049 */
050 @Deprecated
051 protected RandomDataImpl randomData = new RandomDataImpl();
052
053 /**
054 * RNG instance used to generate samples from the distribution.
055 * @since 3.1
056 */
057 protected final RandomGenerator random;
058
059 /** Solver absolute accuracy for inverse cumulative computation */
060 private double solverAbsoluteAccuracy = SOLVER_DEFAULT_ABSOLUTE_ACCURACY;
061
062 /**
063 * @deprecated As of 3.1, to be removed in 4.0. Please use
064 * {@link #AbstractRealDistribution(RandomGenerator)} instead.
065 */
066 @Deprecated
067 protected AbstractRealDistribution() {
068 // Legacy users are only allowed to access the deprecated "randomData".
069 // New users are forbidden to use this constructor.
070 random = null;
071 }
072 /**
073 * @param rng Random number generator.
074 * @since 3.1
075 */
076 protected AbstractRealDistribution(RandomGenerator rng) {
077 random = rng;
078 }
079
080 /**
081 * {@inheritDoc}
082 *
083 * The default implementation uses the identity
084 * <p>{@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}</p>
085 *
086 * @deprecated As of 3.1 (to be removed in 4.0). Please use
087 * {@link #probability(double,double)} instead.
088 */
089 @Deprecated
090 public double cumulativeProbability(double x0, double x1) throws NumberIsTooLargeException {
091 return probability(x0, x1);
092 }
093
094 /**
095 * For a random variable {@code X} whose values are distributed according
096 * to this distribution, this method returns {@code P(x0 < X <= x1)}.
097 *
098 * @param x0 Lower bound (excluded).
099 * @param x1 Upper bound (included).
100 * @return the probability that a random variable with this distribution
101 * takes a value between {@code x0} and {@code x1}, excluding the lower
102 * and including the upper endpoint.
103 * @throws NumberIsTooLargeException if {@code x0 > x1}.
104 *
105 * The default implementation uses the identity
106 * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
107 *
108 * @since 3.1
109 */
110 public double probability(double x0,
111 double x1) {
112 if (x0 > x1) {
113 throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
114 x0, x1, true);
115 }
116 return cumulativeProbability(x1) - cumulativeProbability(x0);
117 }
118
119 /**
120 * {@inheritDoc}
121 *
122 * The default implementation returns
123 * <ul>
124 * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
125 * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
126 * </ul>
127 */
128 public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
129 /*
130 * IMPLEMENTATION NOTES
131 * --------------------
132 * Where applicable, use is made of the one-sided Chebyshev inequality
133 * to bracket the root. This inequality states that
134 * P(X - mu >= k * sig) <= 1 / (1 + k^2),
135 * mu: mean, sig: standard deviation. Equivalently
136 * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
137 * F(mu + k * sig) >= k^2 / (1 + k^2).
138 *
139 * For k = sqrt(p / (1 - p)), we find
140 * F(mu + k * sig) >= p,
141 * and (mu + k * sig) is an upper-bound for the root.
142 *
143 * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
144 * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
145 * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
146 * P(X <= mu - k * sig) <= 1 / (1 + k^2),
147 * F(mu - k * sig) <= 1 / (1 + k^2).
148 *
149 * For k = sqrt((1 - p) / p), we find
150 * F(mu - k * sig) <= p,
151 * and (mu - k * sig) is a lower-bound for the root.
152 *
153 * In cases where the Chebyshev inequality does not apply, geometric
154 * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
155 * the root.
156 */
157 if (p < 0.0 || p > 1.0) {
158 throw new OutOfRangeException(p, 0, 1);
159 }
160
161 double lowerBound = getSupportLowerBound();
162 if (p == 0.0) {
163 return lowerBound;
164 }
165
166 double upperBound = getSupportUpperBound();
167 if (p == 1.0) {
168 return upperBound;
169 }
170
171 final double mu = getNumericalMean();
172 final double sig = FastMath.sqrt(getNumericalVariance());
173 final boolean chebyshevApplies;
174 chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
175 Double.isInfinite(sig) || Double.isNaN(sig));
176
177 if (lowerBound == Double.NEGATIVE_INFINITY) {
178 if (chebyshevApplies) {
179 lowerBound = mu - sig * FastMath.sqrt((1. - p) / p);
180 } else {
181 lowerBound = -1.0;
182 while (cumulativeProbability(lowerBound) >= p) {
183 lowerBound *= 2.0;
184 }
185 }
186 }
187
188 if (upperBound == Double.POSITIVE_INFINITY) {
189 if (chebyshevApplies) {
190 upperBound = mu + sig * FastMath.sqrt(p / (1. - p));
191 } else {
192 upperBound = 1.0;
193 while (cumulativeProbability(upperBound) < p) {
194 upperBound *= 2.0;
195 }
196 }
197 }
198
199 final UnivariateFunction toSolve = new UnivariateFunction() {
200
201 public double value(final double x) {
202 return cumulativeProbability(x) - p;
203 }
204 };
205
206 double x = UnivariateSolverUtils.solve(toSolve,
207 lowerBound,
208 upperBound,
209 getSolverAbsoluteAccuracy());
210
211 if (!isSupportConnected()) {
212 /* Test for plateau. */
213 final double dx = getSolverAbsoluteAccuracy();
214 if (x - dx >= getSupportLowerBound()) {
215 double px = cumulativeProbability(x);
216 if (cumulativeProbability(x - dx) == px) {
217 upperBound = x;
218 while (upperBound - lowerBound > dx) {
219 final double midPoint = 0.5 * (lowerBound + upperBound);
220 if (cumulativeProbability(midPoint) < px) {
221 lowerBound = midPoint;
222 } else {
223 upperBound = midPoint;
224 }
225 }
226 return upperBound;
227 }
228 }
229 }
230 return x;
231 }
232
233 /**
234 * Returns the solver absolute accuracy for inverse cumulative computation.
235 * You can override this method in order to use a Brent solver with an
236 * absolute accuracy different from the default.
237 *
238 * @return the maximum absolute error in inverse cumulative probability estimates
239 */
240 protected double getSolverAbsoluteAccuracy() {
241 return solverAbsoluteAccuracy;
242 }
243
244 /** {@inheritDoc} */
245 public void reseedRandomGenerator(long seed) {
246 random.setSeed(seed);
247 randomData.reSeed(seed);
248 }
249
250 /**
251 * {@inheritDoc}
252 *
253 * The default implementation uses the
254 * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
255 * inversion method.
256 * </a>
257 */
258 public double sample() {
259 return inverseCumulativeProbability(random.nextDouble());
260 }
261
262 /**
263 * {@inheritDoc}
264 *
265 * The default implementation generates the sample by calling
266 * {@link #sample()} in a loop.
267 */
268 public double[] sample(int sampleSize) {
269 if (sampleSize <= 0) {
270 throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SAMPLES,
271 sampleSize);
272 }
273 double[] out = new double[sampleSize];
274 for (int i = 0; i < sampleSize; i++) {
275 out[i] = sample();
276 }
277 return out;
278 }
279
280 /**
281 * {@inheritDoc}
282 *
283 * @return zero.
284 * @since 3.1
285 */
286 public double probability(double x) {
287 return 0d;
288 }
289 }
290