1 /* 2 * Copyright 2003-2004 The Apache Software Foundation. 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); 5 * you may not use this file except in compliance with the License. 6 * You may obtain a copy of the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 * See the License for the specific language governing permissions and 14 * limitations under the License. 15 */ 16 package org.apache.commons.math.distribution; 17 18 import java.io.Serializable; 19 20 import org.apache.commons.math.ConvergenceException; 21 import org.apache.commons.math.FunctionEvaluationException; 22 import org.apache.commons.math.MathException; 23 import org.apache.commons.math.analysis.UnivariateRealFunction; 24 import org.apache.commons.math.analysis.UnivariateRealSolverUtils; 25 26 /** 27 * Base class for continuous distributions. Default implementations are 28 * provided for some of the methods that do not vary from distribution to 29 * distribution. 30 * 31 * @version $Revision: 348519 $ $Date: 2005-11-23 12:12:18 -0700 (Wed, 23 Nov 2005) $ 32 */ 33 public abstract class AbstractContinuousDistribution 34 extends AbstractDistribution 35 implements ContinuousDistribution, Serializable { 36 37 /** Serializable version identifier */ 38 private static final long serialVersionUID = -38038050983108802L; 39 40 /** 41 * Default constructor. 42 */ 43 protected AbstractContinuousDistribution() { 44 super(); 45 } 46 47 /** 48 * For this distribution, X, this method returns the critical point x, such 49 * that P(X < x) = <code>p</code>. 50 * 51 * @param p the desired probability 52 * @return x, such that P(X < x) = <code>p</code> 53 * @throws MathException if the inverse cumulative probability can not be 54 * computed due to convergence or other numerical errors. 55 * @throws IllegalArgumentException if <code>p</code> is not a valid 56 * probability. 57 */ 58 public double inverseCumulativeProbability(final double p) 59 throws MathException { 60 if (p < 0.0 || p > 1.0) { 61 throw new IllegalArgumentException("p must be between 0.0 and 1.0, inclusive."); 62 } 63 64 // by default, do simple root finding using bracketing and default solver. 65 // subclasses can overide if there is a better method. 66 UnivariateRealFunction rootFindingFunction = 67 new UnivariateRealFunction() { 68 69 public double value(double x) throws FunctionEvaluationException { 70 try { 71 return cumulativeProbability(x) - p; 72 } catch (MathException ex) { 73 throw new FunctionEvaluationException 74 (x, "Error computing cdf", ex); 75 } 76 } 77 }; 78 79 // Try to bracket root, test domain endoints if this fails 80 double lowerBound = getDomainLowerBound(p); 81 double upperBound = getDomainUpperBound(p); 82 double[] bracket = null; 83 try { 84 bracket = UnivariateRealSolverUtils.bracket( 85 rootFindingFunction, getInitialDomain(p), 86 lowerBound, upperBound); 87 } catch (ConvergenceException ex) { 88 /* 89 * Check domain endpoints to see if one gives value that is within 90 * the default solver's defaultAbsoluteAccuracy of 0 (will be the 91 * case if density has bounded support and p is 0 or 1). 92 * 93 * TODO: expose the default solver, defaultAbsoluteAccuracy as 94 * a constant. 95 */ 96 if (Math.abs(rootFindingFunction.value(lowerBound)) < 1E-6) { 97 return lowerBound; 98 } 99 if (Math.abs(rootFindingFunction.value(upperBound)) < 1E-6) { 100 return upperBound; 101 } 102 // Failed bracket convergence was not because of corner solution 103 throw new MathException(ex); 104 } 105 106 // find root 107 double root = UnivariateRealSolverUtils.solve(rootFindingFunction, 108 bracket[0],bracket[1]); 109 return root; 110 } 111 112 /** 113 * Access the initial domain value, based on <code>p</code>, used to 114 * bracket a CDF root. This method is used by 115 * {@link #inverseCumulativeProbability(double)} to find critical values. 116 * 117 * @param p the desired probability for the critical value 118 * @return initial domain value 119 */ 120 protected abstract double getInitialDomain(double p); 121 122 /** 123 * Access the domain value lower bound, based on <code>p</code>, used to 124 * bracket a CDF root. This method is used by 125 * {@link #inverseCumulativeProbability(double)} to find critical values. 126 * 127 * @param p the desired probability for the critical value 128 * @return domain value lower bound, i.e. 129 * P(X < <i>lower bound</i>) < <code>p</code> 130 */ 131 protected abstract double getDomainLowerBound(double p); 132 133 /** 134 * Access the domain value upper bound, based on <code>p</code>, used to 135 * bracket a CDF root. This method is used by 136 * {@link #inverseCumulativeProbability(double)} to find critical values. 137 * 138 * @param p the desired probability for the critical value 139 * @return domain value upper bound, i.e. 140 * P(X < <i>upper bound</i>) > <code>p</code> 141 */ 142 protected abstract double getDomainUpperBound(double p); 143 }