View Javadoc

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 &lt; x) = <code>p</code>.
52       *
53       * @param p the desired probability
54       * @return x, such that P(X &lt; 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 &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 }