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: 1.26 $ $Date: 2004/11/07 03:32:48 $
    32   */
    33  public abstract class AbstractContinuousDistribution
    34      extends AbstractDistribution
    35      implements ContinuousDistribution, Serializable {
    36
    37      /*** Serializable version identifier */
    38      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 }