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
019/**
020 * Implementation of the <a href="http://en.wikipedia.org/wiki/Geometric_distribution">geometric distribution</a>.
021 */
022public class GeometricDistribution extends AbstractDiscreteDistribution {
023    /** The probability of success. */
024    private final double probabilityOfSuccess;
025    /** {@code log(p)} where p is the probability of success. */
026    private final double logProbabilityOfSuccess;
027    /** {@code log(1 - p)} where p is the probability of success. */
028    private final double log1mProbabilityOfSuccess;
029
030    /**
031     * Creates a geometric distribution.
032     *
033     * @param p Probability of success.
034     * @throws IllegalArgumentException if {@code p <= 0} or {@code p > 1}.
035     */
036    public GeometricDistribution(double p) {
037        if (p <= 0 || p > 1) {
038            throw new DistributionException(DistributionException.INVALID_PROBABILITY, p);
039        }
040
041        probabilityOfSuccess = p;
042        logProbabilityOfSuccess = Math.log(p);
043        log1mProbabilityOfSuccess = Math.log1p(-p);
044    }
045
046    /**
047     * Access the probability of success for this distribution.
048     *
049     * @return the probability of success.
050     */
051    public double getProbabilityOfSuccess() {
052        return probabilityOfSuccess;
053    }
054
055    /** {@inheritDoc} */
056    @Override
057    public double probability(int x) {
058        if (x < 0) {
059            return 0.0;
060        } else {
061            return Math.exp(log1mProbabilityOfSuccess * x) * probabilityOfSuccess;
062        }
063    }
064
065    /** {@inheritDoc} */
066    @Override
067    public double logProbability(int x) {
068        if (x < 0) {
069            return Double.NEGATIVE_INFINITY;
070        } else {
071            return x * log1mProbabilityOfSuccess + logProbabilityOfSuccess;
072        }
073    }
074
075    /** {@inheritDoc} */
076    @Override
077    public double cumulativeProbability(int x) {
078        if (x < 0) {
079            return 0.0;
080        } else {
081            return -Math.expm1(log1mProbabilityOfSuccess * (x + 1));
082        }
083    }
084
085    /**
086     * {@inheritDoc}
087     *
088     * For probability parameter {@code p}, the mean is {@code (1 - p) / p}.
089     */
090    @Override
091    public double getMean() {
092        return (1 - probabilityOfSuccess) / probabilityOfSuccess;
093    }
094
095    /**
096     * {@inheritDoc}
097     *
098     * For probability parameter {@code p}, the variance is
099     * {@code (1 - p) / (p * p)}.
100     */
101    @Override
102    public double getVariance() {
103        return (1 - probabilityOfSuccess) / (probabilityOfSuccess * probabilityOfSuccess);
104    }
105
106    /**
107     * {@inheritDoc}
108     *
109     * The lower bound of the support is always 0.
110     *
111     * @return lower bound of the support (always 0)
112     */
113    @Override
114    public int getSupportLowerBound() {
115        return 0;
116    }
117
118    /**
119     * {@inheritDoc}
120     *
121     * The upper bound of the support is infinite (which we approximate as
122     * {@code Integer.MAX_VALUE}).
123     *
124     * @return upper bound of the support (always Integer.MAX_VALUE)
125     */
126    @Override
127    public int getSupportUpperBound() {
128        return Integer.MAX_VALUE;
129    }
130
131    /**
132     * {@inheritDoc}
133     *
134     * The support of this distribution is connected.
135     *
136     * @return {@code true}
137     */
138    @Override
139    public boolean isSupportConnected() {
140        return true;
141    }
142
143    /**
144     * {@inheritDoc}
145     */
146    @Override
147    public int inverseCumulativeProbability(double p) {
148        if (p < 0 ||
149            p > 1) {
150            throw new DistributionException(DistributionException.INVALID_PROBABILITY, p);
151        }
152        if (p == 1) {
153            return Integer.MAX_VALUE;
154        }
155        if (p == 0) {
156            return 0;
157        }
158        return Math.max(0, (int) Math.ceil(Math.log1p(-p) / log1mProbabilityOfSuccess - 1));
159    }
160}