001: /*
002: * Copyright 2003-2004 The Apache Software Foundation.
003: *
004: * Licensed under the Apache License, Version 2.0 (the "License");
005: * you may not use this file except in compliance with the License.
006: * You may obtain a copy of the License at
007: *
008: * http://www.apache.org/licenses/LICENSE-2.0
009: *
010: * Unless required by applicable law or agreed to in writing, software
011: * distributed under the License is distributed on an "AS IS" BASIS,
012: * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
013: * See the License for the specific language governing permissions and
014: * limitations under the License.
015: */
016: package org.apache.commons.math.distribution;
017:
018: import java.io.Serializable;
019:
020: import org.apache.commons.math.ConvergenceException;
021: import org.apache.commons.math.FunctionEvaluationException;
022: import org.apache.commons.math.MathException;
023: import org.apache.commons.math.analysis.UnivariateRealFunction;
024: import org.apache.commons.math.analysis.UnivariateRealSolverUtils;
025:
026: /**
027: * Base class for continuous distributions. Default implementations are
028: * provided for some of the methods that do not vary from distribution to
029: * distribution.
030: *
031: * @version $Revision: 348519 $ $Date: 2005-11-23 12:12:18 -0700 (Wed, 23 Nov 2005) $
032: */
033: public abstract class AbstractContinuousDistribution extends
034: AbstractDistribution implements ContinuousDistribution,
035: Serializable {
036:
037: /** Serializable version identifier */
038: private static final long serialVersionUID = -38038050983108802L;
039:
040: /**
041: * Default constructor.
042: */
043: protected AbstractContinuousDistribution() {
044: super ();
045: }
046:
047: /**
048: * For this distribution, X, this method returns the critical point x, such
049: * that P(X < x) = <code>p</code>.
050: *
051: * @param p the desired probability
052: * @return x, such that P(X < x) = <code>p</code>
053: * @throws MathException if the inverse cumulative probability can not be
054: * computed due to convergence or other numerical errors.
055: * @throws IllegalArgumentException if <code>p</code> is not a valid
056: * probability.
057: */
058: public double inverseCumulativeProbability(final double p)
059: throws MathException {
060: if (p < 0.0 || p > 1.0) {
061: throw new IllegalArgumentException(
062: "p must be between 0.0 and 1.0, inclusive.");
063: }
064:
065: // by default, do simple root finding using bracketing and default solver.
066: // subclasses can overide if there is a better method.
067: UnivariateRealFunction rootFindingFunction = new UnivariateRealFunction() {
068:
069: public double value(double x)
070: throws FunctionEvaluationException {
071: try {
072: return cumulativeProbability(x) - p;
073: } catch (MathException ex) {
074: throw new FunctionEvaluationException(x,
075: "Error computing cdf", ex);
076: }
077: }
078: };
079:
080: // Try to bracket root, test domain endoints if this fails
081: double lowerBound = getDomainLowerBound(p);
082: double upperBound = getDomainUpperBound(p);
083: double[] bracket = null;
084: try {
085: bracket = UnivariateRealSolverUtils.bracket(
086: rootFindingFunction, getInitialDomain(p),
087: lowerBound, upperBound);
088: } catch (ConvergenceException ex) {
089: /*
090: * Check domain endpoints to see if one gives value that is within
091: * the default solver's defaultAbsoluteAccuracy of 0 (will be the
092: * case if density has bounded support and p is 0 or 1).
093: *
094: * TODO: expose the default solver, defaultAbsoluteAccuracy as
095: * a constant.
096: */
097: if (Math.abs(rootFindingFunction.value(lowerBound)) < 1E-6) {
098: return lowerBound;
099: }
100: if (Math.abs(rootFindingFunction.value(upperBound)) < 1E-6) {
101: return upperBound;
102: }
103: // Failed bracket convergence was not because of corner solution
104: throw new MathException(ex);
105: }
106:
107: // find root
108: double root = UnivariateRealSolverUtils.solve(
109: rootFindingFunction, bracket[0], bracket[1]);
110: return root;
111: }
112:
113: /**
114: * Access the initial domain value, based on <code>p</code>, used to
115: * bracket a CDF root. This method is used by
116: * {@link #inverseCumulativeProbability(double)} to find critical values.
117: *
118: * @param p the desired probability for the critical value
119: * @return initial domain value
120: */
121: protected abstract double getInitialDomain(double p);
122:
123: /**
124: * Access the domain value lower bound, based on <code>p</code>, used to
125: * bracket a CDF root. This method is used by
126: * {@link #inverseCumulativeProbability(double)} to find critical values.
127: *
128: * @param p the desired probability for the critical value
129: * @return domain value lower bound, i.e.
130: * P(X < <i>lower bound</i>) < <code>p</code>
131: */
132: protected abstract double getDomainLowerBound(double p);
133:
134: /**
135: * Access the domain value upper bound, based on <code>p</code>, used to
136: * bracket a CDF root. This method is used by
137: * {@link #inverseCumulativeProbability(double)} to find critical values.
138: *
139: * @param p the desired probability for the critical value
140: * @return domain value upper bound, i.e.
141: * P(X < <i>upper bound</i>) > <code>p</code>
142: */
143: protected abstract double getDomainUpperBound(double p);
144: }
|