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 */
017package org.apache.commons.statistics.distribution;
018
019import org.apache.commons.numbers.gamma.RegularizedGamma;
020import org.apache.commons.rng.UniformRandomProvider;
021import org.apache.commons.rng.sampling.distribution.PoissonSampler;
022
023/**
024 * Implementation of the <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution</a>.
025 */
026public class PoissonDistribution extends AbstractDiscreteDistribution {
027    /** ln(2 &pi;). */
028    private static final double LOG_TWO_PI = Math.log(2 * Math.PI);
029    /** Default maximum number of iterations. */
030    private static final int DEFAULT_MAX_ITERATIONS = 10000000;
031    /** Default convergence criterion. */
032    private static final double DEFAULT_EPSILON = 1e-12;
033    /** Distribution used to compute normal approximation. */
034    private final NormalDistribution normal;
035    /** Mean of the distribution. */
036    private final double mean;
037    /** Maximum number of iterations for cumulative probability. */
038    private final int maxIterations;
039    /** Convergence criterion for cumulative probability. */
040    private final double epsilon;
041
042    /**
043     * Creates a new Poisson distribution with specified mean.
044     *
045     * @param p the Poisson mean
046     * @throws IllegalArgumentException if {@code p <= 0}.
047     */
048    public PoissonDistribution(double p) {
049        this(p, DEFAULT_EPSILON, DEFAULT_MAX_ITERATIONS);
050    }
051
052    /**
053     * Creates a new Poisson distribution with specified mean, convergence
054     * criterion and maximum number of iterations.
055     *
056     * @param p Poisson mean.
057     * @param epsilon Convergence criterion for cumulative probabilities.
058     * @param maxIterations Maximum number of iterations for cumulative
059     * probabilities.
060     * @throws IllegalArgumentException if {@code p <= 0}.
061     */
062    private PoissonDistribution(double p,
063                                double epsilon,
064                                int maxIterations) {
065        if (p <= 0) {
066            throw new DistributionException(DistributionException.NEGATIVE, p);
067        }
068        mean = p;
069        this.epsilon = epsilon;
070        this.maxIterations = maxIterations;
071
072        normal = new NormalDistribution(p, Math.sqrt(p));
073    }
074
075    /** {@inheritDoc} */
076    @Override
077    public double probability(int x) {
078        final double logProbability = logProbability(x);
079        return logProbability == Double.NEGATIVE_INFINITY ? 0 : Math.exp(logProbability);
080    }
081
082    /** {@inheritDoc} */
083    @Override
084    public double logProbability(int x) {
085        double ret;
086        if (x < 0 || x == Integer.MAX_VALUE) {
087            ret = Double.NEGATIVE_INFINITY;
088        } else if (x == 0) {
089            ret = -mean;
090        } else {
091            ret = -SaddlePointExpansionUtils.getStirlingError(x) -
092                  SaddlePointExpansionUtils.getDeviancePart(x, mean) -
093                  0.5 * LOG_TWO_PI - 0.5 * Math.log(x);
094        }
095        return ret;
096    }
097
098    /** {@inheritDoc} */
099    @Override
100    public double cumulativeProbability(int x) {
101        if (x < 0) {
102            return 0;
103        }
104        if (x == Integer.MAX_VALUE) {
105            return 1;
106        }
107        return RegularizedGamma.Q.value((double) x + 1, mean, epsilon,
108                                        maxIterations);
109    }
110
111    /**
112     * Calculates the Poisson distribution function using a normal
113     * approximation. The {@code N(mean, sqrt(mean))} distribution is used
114     * to approximate the Poisson distribution. The computation uses
115     * "half-correction" (evaluating the normal distribution function at
116     * {@code x + 0.5}).
117     *
118     * @param x Upper bound, inclusive.
119     * @return the distribution function value calculated using a normal
120     * approximation.
121     */
122    public double normalApproximateProbability(int x)  {
123        // Calculate the probability using half-correction.
124        return normal.cumulativeProbability(x + 0.5);
125    }
126
127    /** {@inheritDoc} */
128    @Override
129    public double getMean() {
130        return mean;
131    }
132
133    /**
134     * {@inheritDoc}
135     *
136     * The variance is equal to the {@link #getMean() mean}.
137     */
138    @Override
139    public double getVariance() {
140        return getMean();
141    }
142
143    /**
144     * {@inheritDoc}
145     *
146     * The lower bound of the support is always 0 no matter the mean parameter.
147     *
148     * @return lower bound of the support (always 0)
149     */
150    @Override
151    public int getSupportLowerBound() {
152        return 0;
153    }
154
155    /**
156     * {@inheritDoc}
157     *
158     * The upper bound of the support is positive infinity,
159     * regardless of the parameter values. There is no integer infinity,
160     * so this method returns {@code Integer.MAX_VALUE}.
161     *
162     * @return upper bound of the support (always {@code Integer.MAX_VALUE} for
163     * positive infinity)
164     */
165    @Override
166    public int getSupportUpperBound() {
167        return Integer.MAX_VALUE;
168    }
169
170    /**
171     * {@inheritDoc}
172     *
173     * The support of this distribution is connected.
174     *
175     * @return {@code true}
176     */
177    @Override
178    public boolean isSupportConnected() {
179        return true;
180    }
181
182    /**{@inheritDoc} */
183    @Override
184    public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
185        // Poisson distribution sampler.
186        return new PoissonSampler(rng, mean)::sample;
187    }
188}