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 &lt; x) = <code>p</code>.
071         *
072         * @param p the desired probability
073         * @return x, such that P(X &lt; 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 &lt; <i>lower bound</i>) &lt; <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 &lt; <i>upper bound</i>) &gt; <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    }