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.ConvergenceException;
022 import org.apache.commons.math.MathException;
023 import org.apache.commons.math.MathRuntimeException;
024 import org.apache.commons.math.analysis.UnivariateRealFunction;
025 import org.apache.commons.math.analysis.solvers.BrentSolver;
026 import org.apache.commons.math.analysis.solvers.UnivariateRealSolverUtils;
027 import org.apache.commons.math.FunctionEvaluationException;
028 import org.apache.commons.math.exception.util.LocalizedFormats;
029 import org.apache.commons.math.random.RandomDataImpl;
030 import org.apache.commons.math.util.FastMath;
031
032 /**
033 * Base class for continuous distributions. Default implementations are
034 * provided for some of the methods that do not vary from distribution to
035 * distribution.
036 *
037 * @version $Revision: 1073498 $ $Date: 2011-02-22 21:57:26 +0100 (mar. 22 f??vr. 2011) $
038 */
039 public abstract class AbstractContinuousDistribution
040 extends AbstractDistribution
041 implements ContinuousDistribution, Serializable {
042
043 /** Serializable version identifier */
044 private static final long serialVersionUID = -38038050983108802L;
045
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 /**
053 * Solver absolute accuracy for inverse cumulative computation
054 * @since 2.1
055 */
056 private double solverAbsoluteAccuracy = BrentSolver.DEFAULT_ABSOLUTE_ACCURACY;
057
058 /**
059 * Default constructor.
060 */
061 protected AbstractContinuousDistribution() {
062 super();
063 }
064
065 /**
066 * Return the probability density for a particular point.
067 * @param x The point at which the density should be computed.
068 * @return The pdf at point x.
069 * @throws MathRuntimeException if the specialized class hasn't implemented this function
070 * @since 2.1
071 */
072 public double density(double x) throws MathRuntimeException {
073 throw new MathRuntimeException(new UnsupportedOperationException(),
074 LocalizedFormats.NO_DENSITY_FOR_THIS_DISTRIBUTION);
075 }
076
077 /**
078 * For this distribution, X, this method returns the critical point x, such
079 * that P(X < x) = <code>p</code>.
080 *
081 * @param p the desired probability
082 * @return x, such that P(X < x) = <code>p</code>
083 * @throws MathException if the inverse cumulative probability can not be
084 * computed due to convergence or other numerical errors.
085 * @throws IllegalArgumentException if <code>p</code> is not a valid
086 * probability.
087 */
088 public double inverseCumulativeProbability(final double p)
089 throws MathException {
090 if (p < 0.0 || p > 1.0) {
091 throw MathRuntimeException.createIllegalArgumentException(
092 LocalizedFormats.OUT_OF_RANGE_SIMPLE, p, 0.0, 1.0);
093 }
094
095 // by default, do simple root finding using bracketing and default solver.
096 // subclasses can override if there is a better method.
097 UnivariateRealFunction rootFindingFunction =
098 new UnivariateRealFunction() {
099 public double value(double x) throws FunctionEvaluationException {
100 double ret = Double.NaN;
101 try {
102 ret = cumulativeProbability(x) - p;
103 } catch (MathException ex) {
104 throw new FunctionEvaluationException(x, ex.getSpecificPattern(), ex.getGeneralPattern(), ex.getArguments());
105 }
106 if (Double.isNaN(ret)) {
107 throw new FunctionEvaluationException(x, LocalizedFormats.CUMULATIVE_PROBABILITY_RETURNED_NAN, x, p);
108 }
109 return ret;
110 }
111 };
112
113 // Try to bracket root, test domain endpoints if this fails
114 double lowerBound = getDomainLowerBound(p);
115 double upperBound = getDomainUpperBound(p);
116 double[] bracket = null;
117 try {
118 bracket = UnivariateRealSolverUtils.bracket(
119 rootFindingFunction, getInitialDomain(p),
120 lowerBound, upperBound);
121 } catch (ConvergenceException ex) {
122 /*
123 * Check domain endpoints to see if one gives value that is within
124 * the default solver's defaultAbsoluteAccuracy of 0 (will be the
125 * case if density has bounded support and p is 0 or 1).
126 */
127 if (FastMath.abs(rootFindingFunction.value(lowerBound)) < getSolverAbsoluteAccuracy()) {
128 return lowerBound;
129 }
130 if (FastMath.abs(rootFindingFunction.value(upperBound)) < getSolverAbsoluteAccuracy()) {
131 return upperBound;
132 }
133 // Failed bracket convergence was not because of corner solution
134 throw new MathException(ex);
135 }
136
137 // find root
138 double root = UnivariateRealSolverUtils.solve(rootFindingFunction,
139 // override getSolverAbsoluteAccuracy() to use a Brent solver with
140 // absolute accuracy different from BrentSolver default
141 bracket[0],bracket[1], getSolverAbsoluteAccuracy());
142 return root;
143 }
144
145 /**
146 * Reseeds the random generator used to generate samples.
147 *
148 * @param seed the new seed
149 * @since 2.2
150 */
151 public void reseedRandomGenerator(long seed) {
152 randomData.reSeed(seed);
153 }
154
155 /**
156 * Generates a random value sampled from this distribution. The default
157 * implementation uses the
158 * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling"> inversion method.</a>
159 *
160 * @return random value
161 * @since 2.2
162 * @throws MathException if an error occurs generating the random value
163 */
164 public double sample() throws MathException {
165 return randomData.nextInversionDeviate(this);
166 }
167
168 /**
169 * Generates a random sample from the distribution. The default implementation
170 * generates the sample by calling {@link #sample()} in a loop.
171 *
172 * @param sampleSize number of random values to generate
173 * @since 2.2
174 * @return an array representing the random sample
175 * @throws MathException if an error occurs generating the sample
176 * @throws IllegalArgumentException if sampleSize is not positive
177 */
178 public double[] sample(int sampleSize) throws MathException {
179 if (sampleSize <= 0) {
180 MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NOT_POSITIVE_SAMPLE_SIZE, sampleSize);
181 }
182 double[] out = new double[sampleSize];
183 for (int i = 0; i < sampleSize; i++) {
184 out[i] = sample();
185 }
186 return out;
187 }
188
189 /**
190 * Access the initial domain value, based on <code>p</code>, used to
191 * bracket a CDF root. This method is used by
192 * {@link #inverseCumulativeProbability(double)} to find critical values.
193 *
194 * @param p the desired probability for the critical value
195 * @return initial domain value
196 */
197 protected abstract double getInitialDomain(double p);
198
199 /**
200 * Access the domain value lower bound, based on <code>p</code>, used to
201 * bracket a CDF root. This method is used by
202 * {@link #inverseCumulativeProbability(double)} to find critical values.
203 *
204 * @param p the desired probability for the critical value
205 * @return domain value lower bound, i.e.
206 * P(X < <i>lower bound</i>) < <code>p</code>
207 */
208 protected abstract double getDomainLowerBound(double p);
209
210 /**
211 * Access the domain value upper bound, based on <code>p</code>, used to
212 * bracket a CDF root. This method is used by
213 * {@link #inverseCumulativeProbability(double)} to find critical values.
214 *
215 * @param p the desired probability for the critical value
216 * @return domain value upper bound, i.e.
217 * P(X < <i>upper bound</i>) > <code>p</code>
218 */
219 protected abstract double getDomainUpperBound(double p);
220
221 /**
222 * Returns the solver absolute accuracy for inverse cumulative computation.
223 *
224 * @return the maximum absolute error in inverse cumulative probability estimates
225 * @since 2.1
226 */
227 protected double getSolverAbsoluteAccuracy() {
228 return solverAbsoluteAccuracy;
229 }
230
231 }