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 }