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.UnivariateRealSolverUtils;
027
028 /**
029 * Base class for continuous distributions. Default implementations are
030 * provided for some of the methods that do not vary from distribution to
031 * distribution.
032 *
033 * @version $Revision: 791766 $ $Date: 2009-07-07 05:19:46 -0400 (Tue, 07 Jul 2009) $
034 */
035 public abstract class AbstractContinuousDistribution
036 extends AbstractDistribution
037 implements ContinuousDistribution, Serializable {
038
039 /** Serializable version identifier */
040 private static final long serialVersionUID = -38038050983108802L;
041
042 /**
043 * Default constructor.
044 */
045 protected AbstractContinuousDistribution() {
046 super();
047 }
048
049 /**
050 * For this distribution, X, this method returns the critical point x, such
051 * that P(X < x) = <code>p</code>.
052 *
053 * @param p the desired probability
054 * @return x, such that P(X < x) = <code>p</code>
055 * @throws MathException if the inverse cumulative probability can not be
056 * computed due to convergence or other numerical errors.
057 * @throws IllegalArgumentException if <code>p</code> is not a valid
058 * probability.
059 */
060 public double inverseCumulativeProbability(final double p)
061 throws MathException {
062 if (p < 0.0 || p > 1.0) {
063 throw MathRuntimeException.createIllegalArgumentException(
064 "{0} out of [{1}, {2}] range", p, 0.0, 1.0);
065 }
066
067 // by default, do simple root finding using bracketing and default solver.
068 // subclasses can override if there is a better method.
069 UnivariateRealFunction rootFindingFunction =
070 new UnivariateRealFunction() {
071 public double value(double x) throws FunctionEvaluationException {
072 try {
073 return cumulativeProbability(x) - p;
074 } catch (MathException ex) {
075 throw new FunctionEvaluationException(ex, x, ex.getPattern(), ex.getArguments());
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(rootFindingFunction,
109 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 }