1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math.distribution;
18
19 import java.io.Serializable;
20
21 import org.apache.commons.math.ConvergenceException;
22 import org.apache.commons.math.FunctionEvaluationException;
23 import org.apache.commons.math.MathException;
24 import org.apache.commons.math.MathRuntimeException;
25 import org.apache.commons.math.analysis.UnivariateRealFunction;
26 import org.apache.commons.math.analysis.solvers.UnivariateRealSolverUtils;
27
28 /**
29 * Base class for continuous distributions. Default implementations are
30 * provided for some of the methods that do not vary from distribution to
31 * distribution.
32 *
33 * @version $Revision: 791766 $ $Date: 2009-07-07 05:19:46 -0400 (Tue, 07 Jul 2009) $
34 */
35 public abstract class AbstractContinuousDistribution
36 extends AbstractDistribution
37 implements ContinuousDistribution, Serializable {
38
39 /** Serializable version identifier */
40 private static final long serialVersionUID = -38038050983108802L;
41
42 /**
43 * Default constructor.
44 */
45 protected AbstractContinuousDistribution() {
46 super();
47 }
48
49 /**
50 * For this distribution, X, this method returns the critical point x, such
51 * that P(X < x) = <code>p</code>.
52 *
53 * @param p the desired probability
54 * @return x, such that P(X < x) = <code>p</code>
55 * @throws MathException if the inverse cumulative probability can not be
56 * computed due to convergence or other numerical errors.
57 * @throws IllegalArgumentException if <code>p</code> is not a valid
58 * probability.
59 */
60 public double inverseCumulativeProbability(final double p)
61 throws MathException {
62 if (p < 0.0 || p > 1.0) {
63 throw MathRuntimeException.createIllegalArgumentException(
64 "{0} out of [{1}, {2}] range", p, 0.0, 1.0);
65 }
66
67 // by default, do simple root finding using bracketing and default solver.
68 // subclasses can override if there is a better method.
69 UnivariateRealFunction rootFindingFunction =
70 new UnivariateRealFunction() {
71 public double value(double x) throws FunctionEvaluationException {
72 try {
73 return cumulativeProbability(x) - p;
74 } catch (MathException ex) {
75 throw new FunctionEvaluationException(ex, x, ex.getPattern(), ex.getArguments());
76 }
77 }
78 };
79
80 // Try to bracket root, test domain endoints if this fails
81 double lowerBound = getDomainLowerBound(p);
82 double upperBound = getDomainUpperBound(p);
83 double[] bracket = null;
84 try {
85 bracket = UnivariateRealSolverUtils.bracket(
86 rootFindingFunction, getInitialDomain(p),
87 lowerBound, upperBound);
88 } catch (ConvergenceException ex) {
89 /*
90 * Check domain endpoints to see if one gives value that is within
91 * the default solver's defaultAbsoluteAccuracy of 0 (will be the
92 * case if density has bounded support and p is 0 or 1).
93 *
94 * TODO: expose the default solver, defaultAbsoluteAccuracy as
95 * a constant.
96 */
97 if (Math.abs(rootFindingFunction.value(lowerBound)) < 1E-6) {
98 return lowerBound;
99 }
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 }