AbstractIntegerDistribution.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.distribution;

  19. import java.io.Serializable;

  21. import org.apache.commons.math3.exception.MathInternalError;
  22. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  23. import org.apache.commons.math3.exception.NumberIsTooLargeException;
  24. import org.apache.commons.math3.exception.OutOfRangeException;
  25. import org.apache.commons.math3.exception.util.LocalizedFormats;
  26. import org.apache.commons.math3.random.RandomGenerator;
  27. import org.apache.commons.math3.random.RandomDataImpl;
  28. import org.apache.commons.math3.util.FastMath;

  30. /**
  31.  * Base class for integer-valued discrete distributions.  Default
  32.  * implementations are provided for some of the methods that do not vary
  33.  * from distribution to distribution.
  34.  *
  35.  */
  36. public abstract class AbstractIntegerDistribution implements IntegerDistribution, Serializable {

  38.     /** Serializable version identifier */
  39.     private static final long serialVersionUID = -1146319659338487221L;

  41.     /**
  42.      * RandomData instance used to generate samples from the distribution.
  43.      * @deprecated As of 3.1, to be removed in 4.0. Please use the
  44.      * {@link #random} instance variable instead.
  45.      */
  46.     @Deprecated
  47.     protected final RandomDataImpl randomData = new RandomDataImpl();

  49.     /**
  50.      * RNG instance used to generate samples from the distribution.
  51.      * @since 3.1
  52.      */
  53.     protected final RandomGenerator random;

  55.     /**
  56.      * @deprecated As of 3.1, to be removed in 4.0. Please use
  57.      * {@link #AbstractIntegerDistribution(RandomGenerator)} instead.
  58.      */
  59.     @Deprecated
  60.     protected AbstractIntegerDistribution() {
  61.         // Legacy users are only allowed to access the deprecated "randomData".
  62.         // New users are forbidden to use this constructor.
  63.         random = null;
  64.     }

  66.     /**
  67.      * @param rng Random number generator.
  68.      * @since 3.1
  69.      */
  70.     protected AbstractIntegerDistribution(RandomGenerator rng) {
  71.         random = rng;
  72.     }

  74.     /**
  75.      * {@inheritDoc}
  76.      *
  77.      * The default implementation uses the identity
  78.      * <p>{@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}</p>
  79.      */
  80.     public double cumulativeProbability(int x0, int x1) throws NumberIsTooLargeException {
  81.         if (x1 < x0) {
  82.             throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
  83.                     x0, x1, true);
  84.         }
  85.         return cumulativeProbability(x1) - cumulativeProbability(x0);
  86.     }

  88.     /**
  89.      * {@inheritDoc}
  90.      *
  91.      * The default implementation returns
  92.      * <ul>
  93.      * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
  94.      * <li>{@link #getSupportUpperBound()} for {@code p = 1}, and</li>
  95.      * <li>{@link #solveInverseCumulativeProbability(double, int, int)} for
  96.      *     {@code 0 < p < 1}.</li>
  97.      * </ul>
  98.      */
  99.     public int inverseCumulativeProbability(final double p) throws OutOfRangeException {
  100.         if (p < 0.0 || p > 1.0) {
  101.             throw new OutOfRangeException(p, 0, 1);
  102.         }

  104.         int lower = getSupportLowerBound();
  105.         if (p == 0.0) {
  106.             return lower;
  107.         }
  108.         if (lower == Integer.MIN_VALUE) {
  109.             if (checkedCumulativeProbability(lower) >= p) {
  110.                 return lower;
  111.             }
  112.         } else {
  113.             lower -= 1; // this ensures cumulativeProbability(lower) < p, which
  114.                         // is important for the solving step
  115.         }

  117.         int upper = getSupportUpperBound();
  118.         if (p == 1.0) {
  119.             return upper;
  120.         }

  122.         // use the one-sided Chebyshev inequality to narrow the bracket
  123.         // cf. AbstractRealDistribution.inverseCumulativeProbability(double)
  124.         final double mu = getNumericalMean();
  125.         final double sigma = FastMath.sqrt(getNumericalVariance());
  126.         final boolean chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
  127.                 Double.isInfinite(sigma) || Double.isNaN(sigma) || sigma == 0.0);
  128.         if (chebyshevApplies) {
  129.             double k = FastMath.sqrt((1.0 - p) / p);
  130.             double tmp = mu - k * sigma;
  131.             if (tmp > lower) {
  132.                 lower = ((int) FastMath.ceil(tmp)) - 1;
  133.             }
  134.             k = 1.0 / k;
  135.             tmp = mu + k * sigma;
  136.             if (tmp < upper) {
  137.                 upper = ((int) FastMath.ceil(tmp)) - 1;
  138.             }
  139.         }

  141.         return solveInverseCumulativeProbability(p, lower, upper);
  142.     }

  144.     /**
  145.      * This is a utility function used by {@link
  146.      * #inverseCumulativeProbability(double)}. It assumes {@code 0 < p < 1} and
  147.      * that the inverse cumulative probability lies in the bracket {@code
  148.      * (lower, upper]}. The implementation does simple bisection to find the
  149.      * smallest {@code p}-quantile <code>inf{x in Z | P(X<=x) >= p}</code>.
  150.      *
  151.      * @param p the cumulative probability
  152.      * @param lower a value satisfying {@code cumulativeProbability(lower) < p}
  153.      * @param upper a value satisfying {@code p <= cumulativeProbability(upper)}
  154.      * @return the smallest {@code p}-quantile of this distribution
  155.      */
  156.     protected int solveInverseCumulativeProbability(final double p, int lower, int upper) {
  157.         while (lower + 1 < upper) {
  158.             int xm = (lower + upper) / 2;
  159.             if (xm < lower || xm > upper) {
  160.                 /*
  161.                  * Overflow.
  162.                  * There will never be an overflow in both calculation methods
  163.                  * for xm at the same time
  164.                  */
  165.                 xm = lower + (upper - lower) / 2;
  166.             }

  168.             double pm = checkedCumulativeProbability(xm);
  169.             if (pm >= p) {
  170.                 upper = xm;
  171.             } else {
  172.                 lower = xm;
  173.             }
  174.         }
  175.         return upper;
  176.     }

  178.     /** {@inheritDoc} */
  179.     public void reseedRandomGenerator(long seed) {
  180.         random.setSeed(seed);
  181.         randomData.reSeed(seed);
  182.     }

  184.     /**
  185.      * {@inheritDoc}
  186.      *
  187.      * The default implementation uses the
  188.      * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
  189.      * inversion method</a>.
  190.      */
  191.     public int sample() {
  192.         return inverseCumulativeProbability(random.nextDouble());
  193.     }

  195.     /**
  196.      * {@inheritDoc}
  197.      *
  198.      * The default implementation generates the sample by calling
  199.      * {@link #sample()} in a loop.
  200.      */
  201.     public int[] sample(int sampleSize) {
  202.         if (sampleSize <= 0) {
  203.             throw new NotStrictlyPositiveException(
  204.                     LocalizedFormats.NUMBER_OF_SAMPLES, sampleSize);
  205.         }
  206.         int[] out = new int[sampleSize];
  207.         for (int i = 0; i < sampleSize; i++) {
  208.             out[i] = sample();
  209.         }
  210.         return out;
  211.     }

  213.     /**
  214.      * Computes the cumulative probability function and checks for {@code NaN}
  215.      * values returned. Throws {@code MathInternalError} if the value is
  216.      * {@code NaN}. Rethrows any exception encountered evaluating the cumulative
  217.      * probability function. Throws {@code MathInternalError} if the cumulative
  218.      * probability function returns {@code NaN}.
  219.      *
  220.      * @param argument input value
  221.      * @return the cumulative probability
  222.      * @throws MathInternalError if the cumulative probability is {@code NaN}
  223.      */
  224.     private double checkedCumulativeProbability(int argument)
  225.         throws MathInternalError {
  226.         double result = Double.NaN;
  227.         result = cumulativeProbability(argument);
  228.         if (Double.isNaN(result)) {
  229.             throw new MathInternalError(LocalizedFormats
  230.                     .DISCRETE_CUMULATIVE_PROBABILITY_RETURNED_NAN, argument);
  231.         }
  232.         return result;
  233.     }

  235.     /**
  236.      * For a random variable {@code X} whose values are distributed according to
  237.      * this distribution, this method returns {@code log(P(X = x))}, where
  238.      * {@code log} is the natural logarithm. In other words, this method
  239.      * represents the logarithm of the probability mass function (PMF) for the
  240.      * distribution. Note that due to the floating point precision and
  241.      * under/overflow issues, this method will for some distributions be more
  242.      * precise and faster than computing the logarithm of
  243.      * {@link #probability(int)}.
  244.      * <p>
  245.      * The default implementation simply computes the logarithm of {@code probability(x)}.</p>
  246.      *
  247.      * @param x the point at which the PMF is evaluated
  248.      * @return the logarithm of the value of the probability mass function at {@code x}
  249.      */
  250.     public double logProbability(int x) {
  251.         return FastMath.log(probability(x));
  252.     }
  253. }
Created with JaCoCo 0.7.2.201409121644