    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math.distribution;
    018
    019    import java.io.Serializable;
    020
    021    import org.apache.commons.math.MathException;
    022    import org.apache.commons.math.MathRuntimeException;
    023    import org.apache.commons.math.exception.util.LocalizedFormats;
    024    import org.apache.commons.math.special.Beta;
    025    import org.apache.commons.math.util.MathUtils;
    026    import org.apache.commons.math.util.FastMath;
    027
    028    /**
    029     * The default implementation of {@link PascalDistribution}.
    030     * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $
    031     * @since 1.2
    032     */
    033    public class PascalDistributionImpl extends AbstractIntegerDistribution
    034        implements PascalDistribution, Serializable {
    035
    036        /** Serializable version identifier */
    037        private static final long serialVersionUID = 6751309484392813623L;
    038
    039        /** The number of successes */
    040        private int numberOfSuccesses;
    041
    042        /** The probability of success */
    043        private double probabilityOfSuccess;
    044
    045        /**
    046         * Create a Pascal distribution with the given number of trials and
    047         * probability of success.
    048         * @param r the number of successes
    049         * @param p the probability of success
    050         */
    051        public PascalDistributionImpl(int r, double p) {
    052            super();
    053            setNumberOfSuccessesInternal(r);
    054            setProbabilityOfSuccessInternal(p);
    055        }
    056
    057        /**
    058         * Access the number of successes for this distribution.
    059         * @return the number of successes
    060         */
    061        public int getNumberOfSuccesses() {
    062            return numberOfSuccesses;
    063        }
    064
    065        /**
    066         * Access the probability of success for this distribution.
    067         * @return the probability of success
    068         */
    069        public double getProbabilityOfSuccess() {
    070            return probabilityOfSuccess;
    071        }
    072
    073        /**
    074         * Change the number of successes for this distribution.
    075         * @param successes the new number of successes
    076         * @throws IllegalArgumentException if <code>successes</code> is not
    077         *         positive.
    078         * @deprecated as of 2.1 (class will become immutable in 3.0)
    079         */
    080        @Deprecated
    081        public void setNumberOfSuccesses(int successes) {
    082            setNumberOfSuccessesInternal(successes);
    083        }
    084
    085        /**
    086         * Change the number of successes for this distribution.
    087         * @param successes the new number of successes
    088         * @throws IllegalArgumentException if <code>successes</code> is not
    089         *         positive.
    090         */
    091        private void setNumberOfSuccessesInternal(int successes) {
    092            if (successes < 0) {
    093                throw MathRuntimeException.createIllegalArgumentException(
    094                      LocalizedFormats.NEGATIVE_NUMBER_OF_SUCCESSES,
    095                      successes);
    096            }
    097            numberOfSuccesses = successes;
    098        }
    099
    100        /**
    101         * Change the probability of success for this distribution.
    102         * @param p the new probability of success
    103         * @throws IllegalArgumentException if <code>p</code> is not a valid
    104         *         probability.
    105         * @deprecated as of 2.1 (class will become immutable in 3.0)
    106         */
    107        @Deprecated
    108        public void setProbabilityOfSuccess(double p) {
    109            setProbabilityOfSuccessInternal(p);
    110        }
    111
    112        /**
    113         * Change the probability of success for this distribution.
    114         * @param p the new probability of success
    115         * @throws IllegalArgumentException if <code>p</code> is not a valid
    116         *         probability.
    117         */
    118        private void setProbabilityOfSuccessInternal(double p) {
    119            if (p < 0.0 || p > 1.0) {
    120                throw MathRuntimeException.createIllegalArgumentException(
    121                      LocalizedFormats.OUT_OF_RANGE_SIMPLE, p, 0.0, 1.0);
    122            }
    123            probabilityOfSuccess = p;
    124        }
    125
    126        /**
    127         * Access the domain value lower bound, based on <code>p</code>, used to
    128         * bracket a PDF root.
    129         * @param p the desired probability for the critical value
    130         * @return domain value lower bound, i.e. P(X &lt; <i>lower bound</i>) &lt;
    131         *         <code>p</code>
    132         */
    133        @Override
    134        protected int getDomainLowerBound(double p) {
    135            return -1;
    136        }
    137
    138        /**
    139         * Access the domain value upper bound, based on <code>p</code>, used to
    140         * bracket a PDF root.
    141         * @param p the desired probability for the critical value
    142         * @return domain value upper bound, i.e. P(X &lt; <i>upper bound</i>) &gt;
    143         *         <code>p</code>
    144         */
    145        @Override
    146        protected int getDomainUpperBound(double p) {
    147            // use MAX - 1 because MAX causes loop
    148            return Integer.MAX_VALUE - 1;
    149        }
    150
    151        /**
    152         * For this distribution, X, this method returns P(X &le; x).
    153         * @param x the value at which the PDF is evaluated
    154         * @return PDF for this distribution
    155         * @throws MathException if the cumulative probability can not be computed
    156         *         due to convergence or other numerical errors
    157         */
    158        @Override
    159        public double cumulativeProbability(int x) throws MathException {
    160            double ret;
    161            if (x < 0) {
    162                ret = 0.0;
    163            } else {
    164                ret = Beta.regularizedBeta(probabilityOfSuccess,
    165                    numberOfSuccesses, x + 1);
    166            }
    167            return ret;
    168        }
    169
    170        /**
    171         * For this distribution, X, this method returns P(X = x).
    172         * @param x the value at which the PMF is evaluated
    173         * @return PMF for this distribution
    174         */
    175        public double probability(int x) {
    176            double ret;
    177            if (x < 0) {
    178                ret = 0.0;
    179            } else {
    180                ret = MathUtils.binomialCoefficientDouble(x +
    181                      numberOfSuccesses - 1, numberOfSuccesses - 1) *
    182                      FastMath.pow(probabilityOfSuccess, numberOfSuccesses) *
    183                      FastMath.pow(1.0 - probabilityOfSuccess, x);
    184            }
    185            return ret;
    186        }
    187
    188        /**
    189         * For this distribution, X, this method returns the largest x, such that
    190         * P(X &le; x) &le; <code>p</code>.
    191         * <p>
    192         * Returns <code>-1</code> for p=0 and <code>Integer.MAX_VALUE</code>
    193         * for p=1.</p>
    194         * @param p the desired probability
    195         * @return the largest x such that P(X &le; x) <= p
    196         * @throws MathException if the inverse cumulative probability can not be
    197         *         computed due to convergence or other numerical errors.
    198         * @throws IllegalArgumentException if p < 0 or p > 1
    199         */
    200        @Override
    201        public int inverseCumulativeProbability(final double p)
    202            throws MathException {
    203            int ret;
    204
    205            // handle extreme values explicitly
    206            if (p == 0) {
    207                ret = -1;
    208            } else if (p == 1) {
    209                ret = Integer.MAX_VALUE;
    210            } else {
    211                ret = super.inverseCumulativeProbability(p);
    212            }
    213
    214            return ret;
    215        }
    216
    217        /**
    218         * Returns the lower bound of the support for the distribution.
    219         *
    220         * The lower bound of the support is always 0 no matter the parameters.
    221         *
    222         * @return lower bound of the support (always 0)
    223         * @since 2.2
    224         */
    225        public int getSupportLowerBound() {
    226            return 0;
    227        }
    228
    229        /**
    230         * Returns the upper bound of the support for the distribution.
    231         *
    232         * The upper bound of the support is always positive infinity
    233         * no matter the parameters. Positive infinity is represented
    234         * by <code>Integer.MAX_VALUE</code> together with
    235         * {@link #isSupportUpperBoundInclusive()} being <code>false</code>
    236         *
    237         * @return upper bound of the support (always <code>Integer.MAX_VALUE</code> for positive infinity)
    238         * @since 2.2
    239         */
    240        public int getSupportUpperBound() {
    241            return Integer.MAX_VALUE;
    242        }
    243
    244        /**
    245         * Returns the mean.
    246         *
    247         * For number of successes <code>r</code> and
    248         * probability of success <code>p</code>, the mean is
    249         * <code>( r * p ) / ( 1 - p )</code>
    250         *
    251         * @return the mean
    252         * @since 2.2
    253         */
    254        public double getNumericalMean() {
    255            final double p = getProbabilityOfSuccess();
    256            final double r = getNumberOfSuccesses();
    257            return ( r * p ) / ( 1 - p );
    258        }
    259
    260        /**
    261         * Returns the variance.
    262         *
    263         * For number of successes <code>r</code> and
    264         * probability of success <code>p</code>, the mean is
    265         * <code>( r * p ) / ( 1 - p )^2</code>
    266         *
    267         * @return the variance
    268         * @since 2.2
    269         */
    270        public double getNumericalVariance() {
    271            final double p = getProbabilityOfSuccess();
    272            final double r = getNumberOfSuccesses();
    273            final double pInv = 1 - p;
    274            return ( r * p ) / (pInv * pInv);
    275        }
    276    }