View Javadoc

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