PoissonDistribution.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.  *      https://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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.distribution.discrete;

  24. import org.hipparchus.distribution.continuous.NormalDistribution;
  25. import org.hipparchus.exception.LocalizedCoreFormats;
  26. import org.hipparchus.exception.MathIllegalArgumentException;
  27. import org.hipparchus.special.Gamma;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.MathUtils;

  31. /**
  32.  * Implementation of the Poisson distribution.
  33.  *
  34.  * @see <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution (Wikipedia)</a>
  35.  * @see <a href="http://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution (MathWorld)</a>
  36.  */
  37. public class PoissonDistribution extends AbstractIntegerDistribution {
  38.     /** Default maximum number of iterations for cumulative probability calculations. */
  39.     public static final int DEFAULT_MAX_ITERATIONS = 10000000;
  40.     /** Default convergence criterion. */
  41.     public static final double DEFAULT_EPSILON = 1e-12;
  42.     /** Serializable version identifier. */
  43.     private static final long serialVersionUID = 20160320L;
  44.     /** Distribution used to compute normal approximation. */
  45.     private final NormalDistribution normal;
  46.     /** Mean of the distribution. */
  47.     private final double mean;

  49.     /**
  50.      * Maximum number of iterations for cumulative probability. Cumulative
  51.      * probabilities are estimated using either Lanczos series approximation
  52.      * of {@link Gamma#regularizedGammaP(double, double, double, int)}
  53.      * or continued fraction approximation of
  54.      * {@link Gamma#regularizedGammaQ(double, double, double, int)}.
  55.      */
  56.     private final int maxIterations;

  58.     /** Convergence criterion for cumulative probability. */
  59.     private final double epsilon;

  61.     /**
  62.      * Creates a new Poisson distribution with specified mean.
  63.      *
  64.      * @param p the Poisson mean
  65.      * @throws MathIllegalArgumentException if {@code p <= 0}.
  66.      */
  67.     public PoissonDistribution(double p) throws MathIllegalArgumentException {
  68.         this(p, DEFAULT_EPSILON, DEFAULT_MAX_ITERATIONS);
  69.     }

  71.     /**
  72.      * Creates a new Poisson distribution with specified mean, convergence
  73.      * criterion and maximum number of iterations.
  74.      *
  75.      * @param p Poisson mean.
  76.      * @param epsilon Convergence criterion for cumulative probabilities.
  77.      * @param maxIterations the maximum number of iterations for cumulative
  78.      * probabilities.
  79.      * @throws MathIllegalArgumentException if {@code p <= 0}.
  80.      */
  81.     public PoissonDistribution(double p, double epsilon, int maxIterations)
  82.         throws MathIllegalArgumentException {
  83.         if (p <= 0) {
  84.             throw new MathIllegalArgumentException(LocalizedCoreFormats.MEAN, p);
  85.         }
  86.         mean = p;
  87.         this.epsilon = epsilon;
  88.         this.maxIterations = maxIterations;

  90.         // Use the same RNG instance as the parent class.
  91.         normal = new NormalDistribution(p, FastMath.sqrt(p));
  92.     }

  94.     /**
  95.      * Creates a new Poisson distribution with the specified mean and
  96.      * convergence criterion.
  97.      *
  98.      * @param p Poisson mean.
  99.      * @param epsilon Convergence criterion for cumulative probabilities.
  100.      * @throws MathIllegalArgumentException if {@code p <= 0}.
  101.      */
  102.     public PoissonDistribution(double p, double epsilon)
  103.         throws MathIllegalArgumentException {
  104.         this(p, epsilon, DEFAULT_MAX_ITERATIONS);
  105.     }

  107.     /**
  108.      * Creates a new Poisson distribution with the specified mean and maximum
  109.      * number of iterations.
  110.      *
  111.      * @param p Poisson mean.
  112.      * @param maxIterations Maximum number of iterations for cumulative probabilities.
  113.      */
  114.     public PoissonDistribution(double p, int maxIterations) {
  115.         this(p, DEFAULT_EPSILON, maxIterations);
  116.     }

  118.     /**
  119.      * Get the mean for the distribution.
  120.      *
  121.      * @return the mean for the distribution.
  122.      */
  123.     public double getMean() {
  124.         return mean;
  125.     }

  127.     /** {@inheritDoc} */
  128.     @Override
  129.     public double probability(int x) {
  130.         final double logProbability = logProbability(x);
  131.         return logProbability == Double.NEGATIVE_INFINITY ? 0 : FastMath.exp(logProbability);
  132.     }

  134.     /** {@inheritDoc} */
  135.     @Override
  136.     public double logProbability(int x) {
  137.         double ret;
  138.         if (x < 0 || x == Integer.MAX_VALUE) {
  139.             ret = Double.NEGATIVE_INFINITY;
  140.         } else if (x == 0) {
  141.             ret = -mean;
  142.         } else {
  143.             ret = -SaddlePointExpansion.getStirlingError(x) -
  144.                   SaddlePointExpansion.getDeviancePart(x, mean) -
  145.                   0.5 * FastMath.log(MathUtils.TWO_PI) - 0.5 * FastMath.log(x);
  146.         }
  147.         return ret;
  148.     }

  150.     /** {@inheritDoc} */
  151.     @Override
  152.     public double cumulativeProbability(int x) {
  153.         if (x < 0) {
  154.             return 0;
  155.         }
  156.         if (x == Integer.MAX_VALUE) {
  157.             return 1;
  158.         }
  159.         return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon,
  160.                                        maxIterations);
  161.     }

  163.     /**
  164.      * Calculates the Poisson distribution function using a normal
  165.      * approximation. The {@code N(mean, sqrt(mean))} distribution is used
  166.      * to approximate the Poisson distribution. The computation uses
  167.      * "half-correction" (evaluating the normal distribution function at
  168.      * {@code x + 0.5}).
  169.      *
  170.      * @param x Upper bound, inclusive.
  171.      * @return the distribution function value calculated using a normal
  172.      * approximation.
  173.      */
  174.     public double normalApproximateProbability(int x)  {
  175.         // calculate the probability using half-correction
  176.         return normal.cumulativeProbability(x + 0.5);
  177.     }

  179.     /**
  180.      * {@inheritDoc}
  181.      *
  182.      * For mean parameter {@code p}, the mean is {@code p}.
  183.      */
  184.     @Override
  185.     public double getNumericalMean() {
  186.         return getMean();
  187.     }

  189.     /**
  190.      * {@inheritDoc}
  191.      *
  192.      * For mean parameter {@code p}, the variance is {@code p}.
  193.      */
  194.     @Override
  195.     public double getNumericalVariance() {
  196.         return getMean();
  197.     }

  199.     /**
  200.      * {@inheritDoc}
  201.      *
  202.      * The lower bound of the support is always 0 no matter the mean parameter.
  203.      *
  204.      * @return lower bound of the support (always 0)
  205.      */
  206.     @Override
  207.     public int getSupportLowerBound() {
  208.         return 0;
  209.     }

  211.     /**
  212.      * {@inheritDoc}
  213.      *
  214.      * The upper bound of the support is positive infinity,
  215.      * regardless of the parameter values. There is no integer infinity,
  216.      * so this method returns {@code Integer.MAX_VALUE}.
  217.      *
  218.      * @return upper bound of the support (always {@code Integer.MAX_VALUE} for
  219.      * positive infinity)
  220.      */
  221.     @Override
  222.     public int getSupportUpperBound() {
  223.         return Integer.MAX_VALUE;
  224.     }

  226.     /**
  227.      * {@inheritDoc}
  228.      *
  229.      * The support of this distribution is connected.
  230.      *
  231.      * @return {@code true}
  232.      */
  233.     @Override
  234.     public boolean isSupportConnected() {
  235.         return true;
  236.     }
  237. }
Created with JaCoCo 0.8.12.202403310830