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.math.distribution;
018
019 import java.io.Serializable;
020
021 import org.apache.commons.math.analysis.UnivariateRealFunction;
022 import org.apache.commons.math.analysis.solvers.UnivariateRealSolverUtils;
023 import org.apache.commons.math.exception.MathInternalError;
024 import org.apache.commons.math.exception.NotStrictlyPositiveException;
025 import org.apache.commons.math.exception.OutOfRangeException;
026 import org.apache.commons.math.exception.util.LocalizedFormats;
027 import org.apache.commons.math.exception.NumberIsTooLargeException;
028 import org.apache.commons.math.random.RandomDataImpl;
029 import org.apache.commons.math.util.FastMath;
030
031 /**
032 * Base class for continuous distributions. Default implementations are
033 * provided for some of the methods that do not vary from distribution to
034 * distribution.
035 *
036 * @version $Id: AbstractContinuousDistribution.java 1178295 2011-10-03 04:36:27Z psteitz $
037 */
038 public abstract class AbstractContinuousDistribution
039 extends AbstractDistribution
040 implements ContinuousDistribution, Serializable {
041
042 /** Default accuracy. */
043 public static final double SOLVER_DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
044 /** Serializable version identifier */
045 private static final long serialVersionUID = -38038050983108802L;
046 /**
047 * RandomData instance used to generate samples from the distribution
048 * @since 2.2
049 */
050 protected final RandomDataImpl randomData = new RandomDataImpl();
051 /**
052 * Solver absolute accuracy for inverse cumulative computation.
053 * @since 2.1
054 */
055 private double solverAbsoluteAccuracy = SOLVER_DEFAULT_ABSOLUTE_ACCURACY;
056 /**
057 * Default constructor.
058 */
059 protected AbstractContinuousDistribution() {}
060
061 /**
062 * {@inheritDoc}
063 */
064 public abstract double density(double x);
065
066 /**
067 * For this distribution, {@code X}, this method returns the critical
068 * point {@code x}, such that {@code P(X < x) = p}.
069 *
070 * @param p Desired probability.
071 * @return {@code x}, such that {@code P(X < x) = p}.
072 * @throws OutOfRangeException if {@code p} is not a valid probability.
073 */
074 public double inverseCumulativeProbability(final double p) {
075
076 if (p < 0.0 || p > 1.0) {
077 throw new OutOfRangeException(p, 0, 1);
078 }
079
080 // by default, do simple root finding using bracketing and default solver.
081 // subclasses can override if there is a better method.
082 UnivariateRealFunction rootFindingFunction =
083 new UnivariateRealFunction() {
084 public double value(double x) {
085 return cumulativeProbability(x) - p;
086 }
087 };
088
089 // Try to bracket root, test domain endpoints if this fails
090 double lowerBound = getDomainLowerBound(p);
091 double upperBound = getDomainUpperBound(p);
092 double[] bracket = null;
093 try {
094 bracket = UnivariateRealSolverUtils.bracket(
095 rootFindingFunction, getInitialDomain(p),
096 lowerBound, upperBound);
097 } catch (NumberIsTooLargeException ex) {
098 /*
099 * Check domain endpoints to see if one gives value that is within
100 * the default solver's defaultAbsoluteAccuracy of 0 (will be the
101 * case if density has bounded support and p is 0 or 1).
102 */
103 if (FastMath.abs(rootFindingFunction.value(lowerBound)) < getSolverAbsoluteAccuracy()) {
104 return lowerBound;
105 }
106 if (FastMath.abs(rootFindingFunction.value(upperBound)) < getSolverAbsoluteAccuracy()) {
107 return upperBound;
108 }
109 // Failed bracket convergence was not because of corner solution
110 throw new MathInternalError(ex);
111 }
112
113 // find root
114 double root = UnivariateRealSolverUtils.solve(rootFindingFunction,
115 // override getSolverAbsoluteAccuracy() to use a Brent solver with
116 // absolute accuracy different from the default.
117 bracket[0],bracket[1], getSolverAbsoluteAccuracy());
118 return root;
119 }
120
121 /**
122 * Reseed the random generator used to generate samples.
123 *
124 * @param seed New seed.
125 * @since 2.2
126 */
127 public void reseedRandomGenerator(long seed) {
128 randomData.reSeed(seed);
129 }
130
131 /**
132 * Generate a random value sampled from this distribution. The default
133 * implementation uses the
134 * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
135 * inversion method.
136 * </a>
137 *
138 * @return a random value.
139 * @since 2.2
140 */
141 public double sample() {
142 return randomData.nextInversionDeviate(this);
143 }
144
145 /**
146 * Generate a random sample from the distribution. The default implementation
147 * generates the sample by calling {@link #sample()} in a loop.
148 *
149 * @param sampleSize Number of random values to generate.
150 * @return an array representing the random sample.
151 * @throws NotStrictlyPositiveException if {@code sampleSize} is not positive.
152 * @since 2.2
153 */
154 public double[] sample(int sampleSize) {
155 if (sampleSize <= 0) {
156 throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SAMPLES,
157 sampleSize);
158 }
159 double[] out = new double[sampleSize];
160 for (int i = 0; i < sampleSize; i++) {
161 out[i] = sample();
162 }
163 return out;
164 }
165
166 /**
167 * Access the initial domain value, based on {@code p}, used to
168 * bracket a CDF root. This method is used by
169 * {@link #inverseCumulativeProbability(double)} to find critical values.
170 *
171 * @param p Desired probability for the critical value.
172 * @return the initial domain value.
173 */
174 protected abstract double getInitialDomain(double p);
175
176 /**
177 * Access the domain value lower bound, based on {@code p}, used to
178 * bracket a CDF root. This method is used by
179 * {@link #inverseCumulativeProbability(double)} to find critical values.
180 *
181 * @param p Desired probability for the critical value.
182 * @return the domain value lower bound, i.e. {@code P(X < 'lower bound') < p}.
183 */
184 protected abstract double getDomainLowerBound(double p);
185
186 /**
187 * Access the domain value upper bound, based on {@code p}, used to
188 * bracket a CDF root. This method is used by
189 * {@link #inverseCumulativeProbability(double)} to find critical values.
190 *
191 * @param p Desired probability for the critical value.
192 * @return the domain value upper bound, i.e. {@code P(X < 'upper bound') > p}.
193 */
194 protected abstract double getDomainUpperBound(double p);
195
196 /**
197 * Returns the solver absolute accuracy for inverse cumulative computation.
198 * You can override this method in order to use a Brent solver with an
199 * absolute accuracy different from the default.
200 *
201 * @return the maximum absolute error in inverse cumulative probability estimates
202 * @since 2.1
203 */
204 protected double getSolverAbsoluteAccuracy() {
205 return solverAbsoluteAccuracy;
206 }
207
208 /**
209 * Access the lower bound of the support.
210 *
211 * @return lower bound of the support (might be Double.NEGATIVE_INFINITY)
212 */
213 public abstract double getSupportLowerBound();
214
215 /**
216 * Access the upper bound of the support.
217 *
218 * @return upper bound of the support (might be Double.POSITIVE_INFINITY)
219 */
220 public abstract double getSupportUpperBound();
221 }