Coverage Report - org.apache.commons.math.distribution.GammaDistributionImpl

Classes in this File Line Coverage Branch Coverage Complexity
GammaDistributionImpl
100% 
100% 
2.1

 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$ $Date: 2005-02-26 05:11:52 -0800 (Sat, 26 Feb 2005) $
 27  
  */
 28  
 public class GammaDistributionImpl extends AbstractContinuousDistribution
 29  
     implements GammaDistribution, Serializable  {
 30  
 
 31  
     /** Serializable version identifier */
 32  
     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  142
         super();
 47  142
         setAlpha(alpha);
 48  134
         setBeta(beta);
 49  130
     }
 50  
     
 51  
     /**
 52  
      * For this disbution, X, this method returns P(X < 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  1580
         if (x <= 0.0) {
 72  2
             ret = 0.0;
 73  
         } else {
 74  1578
             ret = Gamma.regularizedGammaP(getAlpha(), x / getBeta());
 75  
         }
 76  
     
 77  1580
         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  36
         if (p == 0) {
 96  2
             return 0d;
 97  
         }
 98  34
         if (p == 1) {
 99  2
             return Double.POSITIVE_INFINITY;
 100  
         }
 101  32
         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  150
         if (alpha <= 0.0) {
 111  12
             throw new IllegalArgumentException("alpha must be positive");
 112  
         }
 113  138
         this.alpha = alpha;
 114  138
     }
 115  
     
 116  
     /**
 117  
      * Access the shape parameter, alpha
 118  
      * @return alpha.
 119  
      */
 120  
     public double getAlpha() {
 121  1686
         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  138
         if (beta <= 0.0) {
 131  6
             throw new IllegalArgumentException("beta must be positive");
 132  
         }
 133  132
         this.beta = beta;
 134  132
     }
 135  
     
 136  
     /**
 137  
      * Access the scale parameter, beta
 138  
      * @return beta.
 139  
      */
 140  
     public double getBeta() {
 141  1662
         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  28
         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  28
         if (p < .5) {
 175  
             // use mean
 176  12
             ret = getAlpha() * getBeta();
 177  
         } else {
 178  
             // use max value
 179  16
             ret = Double.MAX_VALUE;
 180  
         }
 181  
         
 182  28
         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  28
         if (p < .5) {
 200  
             // use 1/2 mean
 201  12
             ret = getAlpha() * getBeta() * .5;
 202  
         } else {
 203  
             // use mean
 204  16
             ret = getAlpha() * getBeta();
 205  
         }
 206  
         
 207  28
         return ret;
 208  
     }
 209  
 }