    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.random.RandomDataImpl;
    025    import org.apache.commons.math.util.FastMath;
    026
    027
    028    /**
    029     * Base class for integer-valued discrete distributions.  Default
    030     * implementations are provided for some of the methods that do not vary
    031     * from distribution to distribution.
    032     *
    033     * @version $Revision: 1067494 $ $Date: 2011-02-05 20:49:07 +0100 (sam. 05 f??vr. 2011) $
    034     */
    035    public abstract class AbstractIntegerDistribution extends AbstractDistribution
    036        implements IntegerDistribution, Serializable {
    037
    038       /** Serializable version identifier */
    039        private static final long serialVersionUID = -1146319659338487221L;
    040
    041        /**
    042         * RandomData instance used to generate samples from the distribution
    043         * @since 2.2
    044         */
    045        protected final RandomDataImpl randomData = new RandomDataImpl();
    046
    047        /**
    048         * Default constructor.
    049         */
    050        protected AbstractIntegerDistribution() {
    051            super();
    052        }
    053
    054        /**
    055         * For a random variable X whose values are distributed according
    056         * to this distribution, this method returns P(X &le; x).  In other words,
    057         * this method represents the  (cumulative) distribution function, or
    058         * CDF, for this distribution.
    059         * <p>
    060         * If <code>x</code> does not represent an integer value, the CDF is
    061         * evaluated at the greatest integer less than x.
    062         *
    063         * @param x the value at which the distribution function is evaluated.
    064         * @return cumulative probability that a random variable with this
    065         * distribution takes a value less than or equal to <code>x</code>
    066         * @throws MathException if the cumulative probability can not be
    067         * computed due to convergence or other numerical errors.
    068         */
    069        public double cumulativeProbability(double x) throws MathException {
    070            return cumulativeProbability((int) FastMath.floor(x));
    071        }
    072
    073        /**
    074         * For a random variable X whose values are distributed according
    075         * to this distribution, this method returns P(x0 &le; X &le; x1).
    076         *
    077         * @param x0 the (inclusive) lower bound
    078         * @param x1 the (inclusive) upper bound
    079         * @return the probability that a random variable with this distribution
    080         * will take a value between <code>x0</code> and <code>x1</code>,
    081         * including the endpoints.
    082         * @throws MathException if the cumulative probability can not be
    083         * computed due to convergence or other numerical errors.
    084         * @throws IllegalArgumentException if <code>x0 > x1</code>
    085         */
    086        @Override
    087        public double cumulativeProbability(double x0, double x1)
    088            throws MathException {
    089            if (x0 > x1) {
    090                throw MathRuntimeException.createIllegalArgumentException(
    091                      LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT, x0, x1);
    092            }
    093            if (FastMath.floor(x0) < x0) {
    094                return cumulativeProbability(((int) FastMath.floor(x0)) + 1,
    095                   (int) FastMath.floor(x1)); // don't want to count mass below x0
    096            } else { // x0 is mathematical integer, so use as is
    097                return cumulativeProbability((int) FastMath.floor(x0),
    098                    (int) FastMath.floor(x1));
    099            }
    100        }
    101
    102        /**
    103         * For a random variable X whose values are distributed according
    104         * to this distribution, this method returns P(X &le; x).  In other words,
    105         * this method represents the probability distribution function, or PDF,
    106         * for this distribution.
    107         *
    108         * @param x the value at which the PDF is evaluated.
    109         * @return PDF for this distribution.
    110         * @throws MathException if the cumulative probability can not be
    111         *            computed due to convergence or other numerical errors.
    112         */
    113        public abstract double cumulativeProbability(int x) throws MathException;
    114
    115        /**
    116         * For a random variable X whose values are distributed according
    117         * to this distribution, this method returns P(X = x). In other words, this
    118         * method represents the probability mass function,  or PMF, for the distribution.
    119         * <p>
    120         * If <code>x</code> does not represent an integer value, 0 is returned.
    121         *
    122         * @param x the value at which the probability density function is evaluated
    123         * @return the value of the probability density function at x
    124         */
    125        public double probability(double x) {
    126            double fl = FastMath.floor(x);
    127            if (fl == x) {
    128                return this.probability((int) x);
    129            } else {
    130                return 0;
    131            }
    132        }
    133
    134        /**
    135        * For a random variable X whose values are distributed according
    136         * to this distribution, this method returns P(x0 &le; X &le; x1).
    137         *
    138         * @param x0 the inclusive, lower bound
    139         * @param x1 the inclusive, upper bound
    140         * @return the cumulative probability.
    141         * @throws MathException if the cumulative probability can not be
    142         *            computed due to convergence or other numerical errors.
    143         * @throws IllegalArgumentException if x0 > x1
    144         */
    145        public double cumulativeProbability(int x0, int x1) throws MathException {
    146            if (x0 > x1) {
    147                throw MathRuntimeException.createIllegalArgumentException(
    148                      LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT, x0, x1);
    149            }
    150            return cumulativeProbability(x1) - cumulativeProbability(x0 - 1);
    151        }
    152
    153        /**
    154         * For a random variable X whose values are distributed according
    155         * to this distribution, this method returns the largest x, such
    156         * that P(X &le; x) &le; <code>p</code>.
    157         *
    158         * @param p the desired probability
    159         * @return the largest x such that P(X &le; x) <= p
    160         * @throws MathException if the inverse cumulative probability can not be
    161         *            computed due to convergence or other numerical errors.
    162         * @throws IllegalArgumentException if p < 0 or p > 1
    163         */
    164        public int inverseCumulativeProbability(final double p) throws MathException{
    165            if (p < 0.0 || p > 1.0) {
    166                throw MathRuntimeException.createIllegalArgumentException(
    167                      LocalizedFormats.OUT_OF_RANGE_SIMPLE, p, 0.0, 1.0);
    168            }
    169
    170            // by default, do simple bisection.
    171            // subclasses can override if there is a better method.
    172            int x0 = getDomainLowerBound(p);
    173            int x1 = getDomainUpperBound(p);
    174            double pm;
    175            while (x0 < x1) {
    176                int xm = x0 + (x1 - x0) / 2;
    177                pm = checkedCumulativeProbability(xm);
    178                if (pm > p) {
    179                    // update x1
    180                    if (xm == x1) {
    181                        // this can happen with integer division
    182                        // simply decrement x1
    183                        --x1;
    184                    } else {
    185                        // update x1 normally
    186                        x1 = xm;
    187                    }
    188                } else {
    189                    // update x0
    190                    if (xm == x0) {
    191                        // this can happen with integer division
    192                        // simply increment x0
    193                        ++x0;
    194                    } else {
    195                        // update x0 normally
    196                        x0 = xm;
    197                    }
    198                }
    199            }
    200
    201            // insure x0 is the correct critical point
    202            pm = checkedCumulativeProbability(x0);
    203            while (pm > p) {
    204                --x0;
    205                pm = checkedCumulativeProbability(x0);
    206            }
    207
    208            return x0;
    209        }
    210
    211        /**
    212         * Reseeds the random generator used to generate samples.
    213         *
    214         * @param seed the new seed
    215         * @since 2.2
    216         */
    217        public void reseedRandomGenerator(long seed) {
    218            randomData.reSeed(seed);
    219        }
    220
    221        /**
    222         * Generates a random value sampled from this distribution. The default
    223         * implementation uses the
    224         * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling"> inversion method.</a>
    225         *
    226         * @return random value
    227         * @since 2.2
    228         * @throws MathException if an error occurs generating the random value
    229         */
    230        public int sample() throws MathException {
    231            return randomData.nextInversionDeviate(this);
    232        }
    233
    234        /**
    235         * Generates a random sample from the distribution.  The default implementation
    236         * generates the sample by calling {@link #sample()} in a loop.
    237         *
    238         * @param sampleSize number of random values to generate
    239         * @since 2.2
    240         * @return an array representing the random sample
    241         * @throws MathException if an error occurs generating the sample
    242         * @throws IllegalArgumentException if sampleSize is not positive
    243         */
    244        public int[] sample(int sampleSize) throws MathException {
    245            if (sampleSize <= 0) {
    246                MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NOT_POSITIVE_SAMPLE_SIZE, sampleSize);
    247            }
    248            int[] out = new int[sampleSize];
    249            for (int i = 0; i < sampleSize; i++) {
    250                out[i] = sample();
    251            }
    252            return out;
    253        }
    254
    255        /**
    256         * Computes the cumulative probability function and checks for NaN values returned.
    257         * Throws MathException if the value is NaN. Rethrows any MathException encountered
    258         * evaluating the cumulative probability function. Throws
    259         * MathException if the cumulative probability function returns NaN.
    260         *
    261         * @param argument input value
    262         * @return cumulative probability
    263         * @throws MathException if the cumulative probability is NaN
    264         */
    265        private double checkedCumulativeProbability(int argument) throws MathException {
    266            double result = Double.NaN;
    267                result = cumulativeProbability(argument);
    268            if (Double.isNaN(result)) {
    269                throw new MathException(LocalizedFormats.DISCRETE_CUMULATIVE_PROBABILITY_RETURNED_NAN, argument);
    270            }
    271            return result;
    272        }
    273
    274        /**
    275         * Access the domain value lower bound, based on <code>p</code>, used to
    276         * bracket a PDF root.  This method is used by
    277         * {@link #inverseCumulativeProbability(double)} to find critical values.
    278         *
    279         * @param p the desired probability for the critical value
    280         * @return domain value lower bound, i.e.
    281         *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
    282         */
    283        protected abstract int getDomainLowerBound(double p);
    284
    285        /**
    286         * Access the domain value upper bound, based on <code>p</code>, used to
    287         * bracket a PDF root.  This method is used by
    288         * {@link #inverseCumulativeProbability(double)} to find critical values.
    289         *
    290         * @param p the desired probability for the critical value
    291         * @return domain value upper bound, i.e.
    292         *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>
    293         */
    294        protected abstract int getDomainUpperBound(double p);
    295
    296        /**
    297         * Use this method to get information about whether the lower bound
    298         * of the support is inclusive or not. For discrete support,
    299         * only true here is meaningful.
    300         *
    301         * @return true (always but at Integer.MIN_VALUE because of the nature of discrete support)
    302         * @since 2.2
    303         */
    304        public boolean isSupportLowerBoundInclusive() {
    305            return true;
    306        }
    307
    308        /**
    309         * Use this method to get information about whether the upper bound
    310         * of the support is inclusive or not. For discrete support,
    311         * only true here is meaningful.
    312         *
    313         * @return true (always but at Integer.MAX_VALUE because of the nature of discrete support)
    314         * @since 2.2
    315         */
    316        public boolean isSupportUpperBoundInclusive() {
    317            return true;
    318        }
    319    }