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.FunctionEvaluationException;
023 import org.apache.commons.math.MathException;
024 import org.apache.commons.math.MathRuntimeException;
025 import org.apache.commons.math.analysis.UnivariateRealFunction;
026 import org.apache.commons.math.analysis.solvers.BrentSolver;
027 import org.apache.commons.math.analysis.solvers.UnivariateRealSolverUtils;
028
029 /**
030 * Base class for continuous distributions. Default implementations are
031 * provided for some of the methods that do not vary from distribution to
032 * distribution.
033 *
034 * @version $Revision: 925812 $ $Date: 2010-03-21 11:49:31 -0400 (Sun, 21 Mar 2010) $
035 */
036 public abstract class AbstractContinuousDistribution
037 extends AbstractDistribution
038 implements ContinuousDistribution, Serializable {
039
040 /** Serializable version identifier */
041 private static final long serialVersionUID = -38038050983108802L;
042
043 /**
044 * Solver absolute accuracy for inverse cum computation
045 * @since 2.1
046 */
047 private double solverAbsoluteAccuracy = BrentSolver.DEFAULT_ABSOLUTE_ACCURACY;
048
049 /**
050 * Default constructor.
051 */
052 protected AbstractContinuousDistribution() {
053 super();
054 }
055
056 /**
057 * Return the probability density for a particular point.
058 * @param x The point at which the density should be computed.
059 * @return The pdf at point x.
060 * @throws MathRuntimeException if the specialized class hasn't implemented this function
061 * @since 2.1
062 */
063 public double density(double x) throws MathRuntimeException {
064 throw new MathRuntimeException(new UnsupportedOperationException(),
065 "This distribution does not have a density function implemented");
066 }
067
068 /**
069 * For this distribution, X, this method returns the critical point x, such
070 * that P(X < x) = <code>p</code>.
071 *
072 * @param p the desired probability
073 * @return x, such that P(X < x) = <code>p</code>
074 * @throws MathException if the inverse cumulative probability can not be
075 * computed due to convergence or other numerical errors.
076 * @throws IllegalArgumentException if <code>p</code> is not a valid
077 * probability.
078 */
079 public double inverseCumulativeProbability(final double p)
080 throws MathException {
081 if (p < 0.0 || p > 1.0) {
082 throw MathRuntimeException.createIllegalArgumentException(
083 "{0} out of [{1}, {2}] range", p, 0.0, 1.0);
084 }
085
086 // by default, do simple root finding using bracketing and default solver.
087 // subclasses can override if there is a better method.
088 UnivariateRealFunction rootFindingFunction =
089 new UnivariateRealFunction() {
090 public double value(double x) throws FunctionEvaluationException {
091 double ret = Double.NaN;
092 try {
093 ret = cumulativeProbability(x) - p;
094 } catch (MathException ex) {
095 throw new FunctionEvaluationException(ex, x, ex.getPattern(), ex.getArguments());
096 }
097 if (Double.isNaN(ret)) {
098 throw new FunctionEvaluationException(x,
099 "Cumulative probability function returned NaN for argument {0} p = {1}", x, p);
100 }
101 return ret;
102 }
103 };
104
105 // Try to bracket root, test domain endoints if this fails
106 double lowerBound = getDomainLowerBound(p);
107 double upperBound = getDomainUpperBound(p);
108 double[] bracket = null;
109 try {
110 bracket = UnivariateRealSolverUtils.bracket(
111 rootFindingFunction, getInitialDomain(p),
112 lowerBound, upperBound);
113 } catch (ConvergenceException ex) {
114 /*
115 * Check domain endpoints to see if one gives value that is within
116 * the default solver's defaultAbsoluteAccuracy of 0 (will be the
117 * case if density has bounded support and p is 0 or 1).
118 */
119 if (Math.abs(rootFindingFunction.value(lowerBound)) < getSolverAbsoluteAccuracy()) {
120 return lowerBound;
121 }
122 if (Math.abs(rootFindingFunction.value(upperBound)) < getSolverAbsoluteAccuracy()) {
123 return upperBound;
124 }
125 // Failed bracket convergence was not because of corner solution
126 throw new MathException(ex);
127 }
128
129 // find root
130 double root = UnivariateRealSolverUtils.solve(rootFindingFunction,
131 // override getSolverAbsoluteAccuracy() to use a Brent solver with
132 // absolute accuracy different from BrentSolver default
133 bracket[0],bracket[1], getSolverAbsoluteAccuracy());
134 return root;
135 }
136
137 /**
138 * Access the initial domain value, based on <code>p</code>, used to
139 * bracket a CDF root. This method is used by
140 * {@link #inverseCumulativeProbability(double)} to find critical values.
141 *
142 * @param p the desired probability for the critical value
143 * @return initial domain value
144 */
145 protected abstract double getInitialDomain(double p);
146
147 /**
148 * Access the domain value lower bound, based on <code>p</code>, used to
149 * bracket a CDF root. This method is used by
150 * {@link #inverseCumulativeProbability(double)} to find critical values.
151 *
152 * @param p the desired probability for the critical value
153 * @return domain value lower bound, i.e.
154 * P(X < <i>lower bound</i>) < <code>p</code>
155 */
156 protected abstract double getDomainLowerBound(double p);
157
158 /**
159 * Access the domain value upper bound, based on <code>p</code>, used to
160 * bracket a CDF root. This method is used by
161 * {@link #inverseCumulativeProbability(double)} to find critical values.
162 *
163 * @param p the desired probability for the critical value
164 * @return domain value upper bound, i.e.
165 * P(X < <i>upper bound</i>) > <code>p</code>
166 */
167 protected abstract double getDomainUpperBound(double p);
168
169 /**
170 * Returns the solver absolute accuracy for inverse cum computation.
171 *
172 * @return the maximum absolute error in inverse cumulative probability estimates
173 * @since 2.1
174 */
175 protected double getSolverAbsoluteAccuracy() {
176 return solverAbsoluteAccuracy;
177 }
178 }