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.MathException;
    21  import org.apache.commons.math.special.Gamma;
    22
    23  /***
    24   * The default implementation of {@link GammaDistribution}.
    25   *
    26   * @version $Revision: 355770 $ $Date: 2005-12-10 12:48:57 -0700 (Sat, 10 Dec 2005) $
    27   */
    28  public class GammaDistributionImpl extends AbstractContinuousDistribution
    29      implements GammaDistribution, Serializable  {
    30
    31      /*** Serializable version identifier */
    32      private static final long serialVersionUID = -3239549463135430361L;
    33
    34      /*** The shape parameter. */
    35      private double alpha;
    36
    37      /*** The scale parameter. */
    38      private double beta;
    39
    40      /***
    41       * Create a new gamma distribution with the given alpha and beta values.
    42       * @param alpha the shape parameter.
    43       * @param beta the scale parameter.
    44       */
    45      public GammaDistributionImpl(double alpha, double beta) {
    46          super();
    47          setAlpha(alpha);
    48          setBeta(beta);
    49      }
    50
    51      /***
    52       * For this disbution, X, this method returns P(X &lt; x).
    53       *
    54       * The implementation of this method is based on:
    55       * <ul>
    56       * <li>
    57       * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
    58       * Chi-Squared Distribution</a>, equation (9).</li>
    59       * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>.
    60       * Belmont, CA: Duxbury Press.</li>
    61       * </ul>
    62       *
    63       * @param x the value at which the CDF is evaluated.
    64       * @return CDF for this distribution.
    65       * @throws MathException if the cumulative probability can not be
    66       *            computed due to convergence or other numerical errors.
    67       */
    68      public double cumulativeProbability(double x) throws MathException{
    69          double ret;
    70
    71          if (x <= 0.0) {
    72              ret = 0.0;
    73          } else {
    74              ret = Gamma.regularizedGammaP(getAlpha(), x / getBeta());
    75          }
    76
    77          return ret;
    78      }
    79
    80      /***
    81       * For this distribution, X, this method returns the critical point x, such
    82       * that P(X &lt; x) = <code>p</code>.
    83       * <p>
    84       * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.
    85       *
    86       * @param p the desired probability
    87       * @return x, such that P(X &lt; x) = <code>p</code>
    88       * @throws MathException if the inverse cumulative probability can not be
    89       *         computed due to convergence or other numerical errors.
    90       * @throws IllegalArgumentException if <code>p</code> is not a valid
    91       *         probability.
    92       */
    93      public double inverseCumulativeProbability(final double p)
    94      throws MathException {
    95          if (p == 0) {
    96              return 0d;
    97          }
    98          if (p == 1) {
    99              return Double.POSITIVE_INFINITY;
    100         }
    101         return super.inverseCumulativeProbability(p);
    102     }
    103
    104     /***
    105      * Modify the shape parameter, alpha.
    106      * @param alpha the new shape parameter.
    107      * @throws IllegalArgumentException if <code>alpha</code> is not positive.
    108      */
    109     public void setAlpha(double alpha) {
    110         if (alpha <= 0.0) {
    111             throw new IllegalArgumentException("alpha must be positive");
    112         }
    113         this.alpha = alpha;
    114     }
    115
    116     /***
    117      * Access the shape parameter, alpha
    118      * @return alpha.
    119      */
    120     public double getAlpha() {
    121         return alpha;
    122     }
    123
    124     /***
    125      * Modify the scale parameter, beta.
    126      * @param beta the new scale parameter.
    127      * @throws IllegalArgumentException if <code>beta</code> is not positive.
    128      */
    129     public void setBeta(double beta) {
    130         if (beta <= 0.0) {
    131             throw new IllegalArgumentException("beta must be positive");
    132         }
    133         this.beta = beta;
    134     }
    135
    136     /***
    137      * Access the scale parameter, beta
    138      * @return beta.
    139      */
    140     public double getBeta() {
    141         return beta;
    142     }
    143
    144     /***
    145      * Access the domain value lower bound, based on <code>p</code>, used to
    146      * bracket a CDF root.  This method is used by
    147      * {@link #inverseCumulativeProbability(double)} to find critical values.
    148      *
    149      * @param p the desired probability for the critical value
    150      * @return domain value lower bound, i.e.
    151      *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
    152      */
    153     protected double getDomainLowerBound(double p) {
    154         // TODO: try to improve on this estimate
    155         return Double.MIN_VALUE;
    156     }
    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 double getDomainUpperBound(double p) {
    168         // TODO: try to improve on this estimate
    169         // NOTE: gamma is skewed to the left
    170         // NOTE: therefore, P(X < &mu;) > .5
    171
    172         double ret;
    173
    174         if (p < .5) {
    175             // use mean
    176             ret = getAlpha() * getBeta();
    177         } else {
    178             // use max value
    179             ret = Double.MAX_VALUE;
    180         }
    181
    182         return ret;
    183     }
    184
    185     /***
    186      * Access the initial domain value, based on <code>p</code>, used to
    187      * bracket a CDF root.  This method is used by
    188      * {@link #inverseCumulativeProbability(double)} to find critical values.
    189      *
    190      * @param p the desired probability for the critical value
    191      * @return initial domain value
    192      */
    193     protected double getInitialDomain(double p) {
    194         // TODO: try to improve on this estimate
    195         // Gamma is skewed to the left, therefore, P(X < &mu;) > .5
    196
    197         double ret;
    198
    199         if (p < .5) {
    200             // use 1/2 mean
    201             ret = getAlpha() * getBeta() * .5;
    202         } else {
    203             // use mean
    204             ret = getAlpha() * getBeta();
    205         }
    206
    207         return ret;
    208     }
    209 }