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.analysis.UnivariateRealFunction;
25  import org.apache.commons.math.analysis.UnivariateRealSolverUtils;
26  
27  /**
28   * Base class for continuous distributions.  Default implementations are
29   * provided for some of the methods that do not vary from distribution to
30   * distribution.
31   *  
32   * @version $Revision: 506600 $ $Date: 2007-02-12 12:35:59 -0700 (Mon, 12 Feb 2007) $
33   */
34  public abstract class AbstractContinuousDistribution
35      extends AbstractDistribution
36      implements ContinuousDistribution, Serializable {
37  
38      /** Serializable version identifier */
39      private static final long serialVersionUID = -38038050983108802L;
40      
41      /**
42       * Default constructor.
43       */
44      protected AbstractContinuousDistribution() {
45          super();
46      }
47  
48      /**
49       * For this distribution, X, this method returns the critical point x, such
50       * that P(X &lt; x) = <code>p</code>.
51       *
52       * @param p the desired probability
53       * @return x, such that P(X &lt; x) = <code>p</code>
54       * @throws MathException if the inverse cumulative probability can not be
55       *         computed due to convergence or other numerical errors.
56       * @throws IllegalArgumentException if <code>p</code> is not a valid
57       *         probability.
58       */
59      public double inverseCumulativeProbability(final double p)
60          throws MathException {
61          if (p < 0.0 || p > 1.0) {
62              throw new IllegalArgumentException("p must be between 0.0 and 1.0, inclusive.");
63          }
64  
65          // by default, do simple root finding using bracketing and default solver.
66          // subclasses can overide if there is a better method.
67          UnivariateRealFunction rootFindingFunction =
68              new UnivariateRealFunction() {
69  
70              public double value(double x) throws FunctionEvaluationException {
71                  try {
72                      return cumulativeProbability(x) - p;
73                  } catch (MathException ex) {
74                      throw new FunctionEvaluationException(x, ex.getPattern(), ex.getArguments(), 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 }