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.special.Beta;
024 import org.apache.commons.math.util.MathUtils;
025
026 /**
027 * The default implementation of {@link PascalDistribution}.
028 * @version $Revision: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
029 * @since 1.2
030 */
031 public class PascalDistributionImpl extends AbstractIntegerDistribution
032 implements PascalDistribution, Serializable {
033
034 /** Serializable version identifier */
035 private static final long serialVersionUID = 6751309484392813623L;
036
037 /** The number of successes */
038 private int numberOfSuccesses;
039
040 /** The probability of success */
041 private double probabilityOfSuccess;
042
043 /**
044 * Create a binomial distribution with the given number of trials and
045 * probability of success.
046 * @param r the number of successes
047 * @param p the probability of success
048 */
049 public PascalDistributionImpl(int r, double p) {
050 super();
051 setNumberOfSuccesses(r);
052 setProbabilityOfSuccess(p);
053 }
054
055 /**
056 * Access the number of successes for this distribution.
057 * @return the number of successes
058 */
059 public int getNumberOfSuccesses() {
060 return numberOfSuccesses;
061 }
062
063 /**
064 * Access the probability of success for this distribution.
065 * @return the probability of success
066 */
067 public double getProbabilityOfSuccess() {
068 return probabilityOfSuccess;
069 }
070
071 /**
072 * Change the number of successes for this distribution.
073 * @param successes the new number of successes
074 * @throws IllegalArgumentException if <code>successes</code> is not
075 * positive.
076 */
077 public void setNumberOfSuccesses(int successes) {
078 if (successes < 0) {
079 throw MathRuntimeException.createIllegalArgumentException(
080 "number of successes must be non-negative ({0})",
081 successes);
082 }
083 numberOfSuccesses = successes;
084 }
085
086 /**
087 * Change the probability of success for this distribution.
088 * @param p the new probability of success
089 * @throws IllegalArgumentException if <code>p</code> is not a valid
090 * probability.
091 */
092 public void setProbabilityOfSuccess(double p) {
093 if (p < 0.0 || p > 1.0) {
094 throw MathRuntimeException.createIllegalArgumentException(
095 "{0} out of [{1}, {2}] range", p, 0.0, 1.0);
096 }
097 probabilityOfSuccess = p;
098 }
099
100 /**
101 * Access the domain value lower bound, based on <code>p</code>, used to
102 * bracket a PDF root.
103 * @param p the desired probability for the critical value
104 * @return domain value lower bound, i.e. P(X < <i>lower bound</i>) <
105 * <code>p</code>
106 */
107 @Override
108 protected int getDomainLowerBound(double p) {
109 return -1;
110 }
111
112 /**
113 * Access the domain value upper bound, based on <code>p</code>, used to
114 * bracket a PDF root.
115 * @param p the desired probability for the critical value
116 * @return domain value upper bound, i.e. P(X < <i>upper bound</i>) >
117 * <code>p</code>
118 */
119 @Override
120 protected int getDomainUpperBound(double p) {
121 // use MAX - 1 because MAX causes loop
122 return Integer.MAX_VALUE - 1;
123 }
124
125 /**
126 * For this distribution, X, this method returns P(X ≤ x).
127 * @param x the value at which the PDF is evaluated
128 * @return PDF for this distribution
129 * @throws MathException if the cumulative probability can not be computed
130 * due to convergence or other numerical errors
131 */
132 @Override
133 public double cumulativeProbability(int x) throws MathException {
134 double ret;
135 if (x < 0) {
136 ret = 0.0;
137 } else {
138 ret = Beta.regularizedBeta(getProbabilityOfSuccess(),
139 getNumberOfSuccesses(), x + 1);
140 }
141 return ret;
142 }
143
144 /**
145 * For this distribution, X, this method returns P(X = x).
146 * @param x the value at which the PMF is evaluated
147 * @return PMF for this distribution
148 */
149 public double probability(int x) {
150 double ret;
151 if (x < 0) {
152 ret = 0.0;
153 } else {
154 ret = MathUtils.binomialCoefficientDouble(x +
155 getNumberOfSuccesses() - 1, getNumberOfSuccesses() - 1) *
156 Math.pow(getProbabilityOfSuccess(), getNumberOfSuccesses()) *
157 Math.pow(1.0 - getProbabilityOfSuccess(), x);
158 }
159 return ret;
160 }
161
162 /**
163 * For this distribution, X, this method returns the largest x, such that
164 * P(X ≤ x) ≤ <code>p</code>.
165 * <p>
166 * Returns <code>-1</code> for p=0 and <code>Integer.MAX_VALUE</code>
167 * for p=1.</p>
168 * @param p the desired probability
169 * @return the largest x such that P(X ≤ x) <= p
170 * @throws MathException if the inverse cumulative probability can not be
171 * computed due to convergence or other numerical errors.
172 * @throws IllegalArgumentException if p < 0 or p > 1
173 */
174 @Override
175 public int inverseCumulativeProbability(final double p)
176 throws MathException {
177 int ret;
178
179 // handle extreme values explicitly
180 if (p == 0) {
181 ret = -1;
182 } else if (p == 1) {
183 ret = Integer.MAX_VALUE;
184 } else {
185 ret = super.inverseCumulativeProbability(p);
186 }
187
188 return ret;
189 }
190 }