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 &lt; x) = <code>p</code>.
052         *
053         * @param p the desired probability
054         * @return x, such that P(X &lt; 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 &lt; <i>lower bound</i>) &lt; <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 &lt; <i>upper bound</i>) &gt; <code>p</code> 
142         */
143        protected abstract double getDomainUpperBound(double p);
144    }