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.math.distribution; 18 19 import java.io.Serializable; 20 21 import org.apache.commons.math.MathException; 22 import org.apache.commons.math.MathRuntimeException; 23 import org.apache.commons.math.exception.util.LocalizedFormats; 24 import org.apache.commons.math.special.Beta; 25 import org.apache.commons.math.util.MathUtils; 26 import org.apache.commons.math.util.FastMath; 27 28 /** 29 * The default implementation of {@link PascalDistribution}. 30 * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $ 31 * @since 1.2 32 */ 33 public class PascalDistributionImpl extends AbstractIntegerDistribution 34 implements PascalDistribution, Serializable { 35 36 /** Serializable version identifier */ 37 private static final long serialVersionUID = 6751309484392813623L; 38 39 /** The number of successes */ 40 private int numberOfSuccesses; 41 42 /** The probability of success */ 43 private double probabilityOfSuccess; 44 45 /** 46 * Create a Pascal distribution with the given number of trials and 47 * probability of success. 48 * @param r the number of successes 49 * @param p the probability of success 50 */ 51 public PascalDistributionImpl(int r, double p) { 52 super(); 53 setNumberOfSuccessesInternal(r); 54 setProbabilityOfSuccessInternal(p); 55 } 56 57 /** 58 * Access the number of successes for this distribution. 59 * @return the number of successes 60 */ 61 public int getNumberOfSuccesses() { 62 return numberOfSuccesses; 63 } 64 65 /** 66 * Access the probability of success for this distribution. 67 * @return the probability of success 68 */ 69 public double getProbabilityOfSuccess() { 70 return probabilityOfSuccess; 71 } 72 73 /** 74 * Change the number of successes for this distribution. 75 * @param successes the new number of successes 76 * @throws IllegalArgumentException if <code>successes</code> is not 77 * positive. 78 * @deprecated as of 2.1 (class will become immutable in 3.0) 79 */ 80 @Deprecated 81 public void setNumberOfSuccesses(int successes) { 82 setNumberOfSuccessesInternal(successes); 83 } 84 85 /** 86 * Change the number of successes for this distribution. 87 * @param successes the new number of successes 88 * @throws IllegalArgumentException if <code>successes</code> is not 89 * positive. 90 */ 91 private void setNumberOfSuccessesInternal(int successes) { 92 if (successes < 0) { 93 throw MathRuntimeException.createIllegalArgumentException( 94 LocalizedFormats.NEGATIVE_NUMBER_OF_SUCCESSES, 95 successes); 96 } 97 numberOfSuccesses = successes; 98 } 99 100 /** 101 * Change the probability of success for this distribution. 102 * @param p the new probability of success 103 * @throws IllegalArgumentException if <code>p</code> is not a valid 104 * probability. 105 * @deprecated as of 2.1 (class will become immutable in 3.0) 106 */ 107 @Deprecated 108 public void setProbabilityOfSuccess(double p) { 109 setProbabilityOfSuccessInternal(p); 110 } 111 112 /** 113 * Change the probability of success for this distribution. 114 * @param p the new probability of success 115 * @throws IllegalArgumentException if <code>p</code> is not a valid 116 * probability. 117 */ 118 private void setProbabilityOfSuccessInternal(double p) { 119 if (p < 0.0 || p > 1.0) { 120 throw MathRuntimeException.createIllegalArgumentException( 121 LocalizedFormats.OUT_OF_RANGE_SIMPLE, p, 0.0, 1.0); 122 } 123 probabilityOfSuccess = p; 124 } 125 126 /** 127 * Access the domain value lower bound, based on <code>p</code>, used to 128 * bracket a PDF root. 129 * @param p the desired probability for the critical value 130 * @return domain value lower bound, i.e. P(X < <i>lower bound</i>) < 131 * <code>p</code> 132 */ 133 @Override 134 protected int getDomainLowerBound(double p) { 135 return -1; 136 } 137 138 /** 139 * Access the domain value upper bound, based on <code>p</code>, used to 140 * bracket a PDF root. 141 * @param p the desired probability for the critical value 142 * @return domain value upper bound, i.e. P(X < <i>upper bound</i>) > 143 * <code>p</code> 144 */ 145 @Override 146 protected int getDomainUpperBound(double p) { 147 // use MAX - 1 because MAX causes loop 148 return Integer.MAX_VALUE - 1; 149 } 150 151 /** 152 * For this distribution, X, this method returns P(X ≤ x). 153 * @param x the value at which the PDF is evaluated 154 * @return PDF for this distribution 155 * @throws MathException if the cumulative probability can not be computed 156 * due to convergence or other numerical errors 157 */ 158 @Override 159 public double cumulativeProbability(int x) throws MathException { 160 double ret; 161 if (x < 0) { 162 ret = 0.0; 163 } else { 164 ret = Beta.regularizedBeta(probabilityOfSuccess, 165 numberOfSuccesses, x + 1); 166 } 167 return ret; 168 } 169 170 /** 171 * For this distribution, X, this method returns P(X = x). 172 * @param x the value at which the PMF is evaluated 173 * @return PMF for this distribution 174 */ 175 public double probability(int x) { 176 double ret; 177 if (x < 0) { 178 ret = 0.0; 179 } else { 180 ret = MathUtils.binomialCoefficientDouble(x + 181 numberOfSuccesses - 1, numberOfSuccesses - 1) * 182 FastMath.pow(probabilityOfSuccess, numberOfSuccesses) * 183 FastMath.pow(1.0 - probabilityOfSuccess, x); 184 } 185 return ret; 186 } 187 188 /** 189 * For this distribution, X, this method returns the largest x, such that 190 * P(X ≤ x) ≤ <code>p</code>. 191 * <p> 192 * Returns <code>-1</code> for p=0 and <code>Integer.MAX_VALUE</code> 193 * for p=1.</p> 194 * @param p the desired probability 195 * @return the largest x such that P(X ≤ x) <= p 196 * @throws MathException if the inverse cumulative probability can not be 197 * computed due to convergence or other numerical errors. 198 * @throws IllegalArgumentException if p < 0 or p > 1 199 */ 200 @Override 201 public int inverseCumulativeProbability(final double p) 202 throws MathException { 203 int ret; 204 205 // handle extreme values explicitly 206 if (p == 0) { 207 ret = -1; 208 } else if (p == 1) { 209 ret = Integer.MAX_VALUE; 210 } else { 211 ret = super.inverseCumulativeProbability(p); 212 } 213 214 return ret; 215 } 216 217 /** 218 * Returns the lower bound of the support for the distribution. 219 * 220 * The lower bound of the support is always 0 no matter the parameters. 221 * 222 * @return lower bound of the support (always 0) 223 * @since 2.2 224 */ 225 public int getSupportLowerBound() { 226 return 0; 227 } 228 229 /** 230 * Returns the upper bound of the support for the distribution. 231 * 232 * The upper bound of the support is always positive infinity 233 * no matter the parameters. Positive infinity is represented 234 * by <code>Integer.MAX_VALUE</code> together with 235 * {@link #isSupportUpperBoundInclusive()} being <code>false</code> 236 * 237 * @return upper bound of the support (always <code>Integer.MAX_VALUE</code> for positive infinity) 238 * @since 2.2 239 */ 240 public int getSupportUpperBound() { 241 return Integer.MAX_VALUE; 242 } 243 244 /** 245 * Returns the mean. 246 * 247 * For number of successes <code>r</code> and 248 * probability of success <code>p</code>, the mean is 249 * <code>( r * p ) / ( 1 - p )</code> 250 * 251 * @return the mean 252 * @since 2.2 253 */ 254 public double getNumericalMean() { 255 final double p = getProbabilityOfSuccess(); 256 final double r = getNumberOfSuccesses(); 257 return ( r * p ) / ( 1 - p ); 258 } 259 260 /** 261 * Returns the variance. 262 * 263 * For number of successes <code>r</code> and 264 * probability of success <code>p</code>, the mean is 265 * <code>( r * p ) / ( 1 - p )^2</code> 266 * 267 * @return the variance 268 * @since 2.2 269 */ 270 public double getNumericalVariance() { 271 final double p = getProbabilityOfSuccess(); 272 final double r = getNumberOfSuccesses(); 273 final double pInv = 1 - p; 274 return ( r * p ) / (pInv * pInv); 275 } 276 } 277