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.random.RandomDataImpl; 25 import org.apache.commons.math.util.FastMath; 26 27 28 /** 29 * Base class for integer-valued discrete distributions. Default 30 * implementations are provided for some of the methods that do not vary 31 * from distribution to distribution. 32 * 33 * @version $Revision: 1067494 $ $Date: 2011-02-05 20:49:07 +0100 (sam. 05 fvr. 2011) $ 34 */ 35 public abstract class AbstractIntegerDistribution extends AbstractDistribution 36 implements IntegerDistribution, Serializable { 37 38 /** Serializable version identifier */ 39 private static final long serialVersionUID = -1146319659338487221L; 40 41 /** 42 * RandomData instance used to generate samples from the distribution 43 * @since 2.2 44 */ 45 protected final RandomDataImpl randomData = new RandomDataImpl(); 46 47 /** 48 * Default constructor. 49 */ 50 protected AbstractIntegerDistribution() { 51 super(); 52 } 53 54 /** 55 * For a random variable X whose values are distributed according 56 * to this distribution, this method returns P(X ≤ x). In other words, 57 * this method represents the (cumulative) distribution function, or 58 * CDF, for this distribution. 59 * <p> 60 * If <code>x</code> does not represent an integer value, the CDF is 61 * evaluated at the greatest integer less than x. 62 * 63 * @param x the value at which the distribution function is evaluated. 64 * @return cumulative probability that a random variable with this 65 * distribution takes a value less than or equal to <code>x</code> 66 * @throws MathException if the cumulative probability can not be 67 * computed due to convergence or other numerical errors. 68 */ 69 public double cumulativeProbability(double x) throws MathException { 70 return cumulativeProbability((int) FastMath.floor(x)); 71 } 72 73 /** 74 * For a random variable X whose values are distributed according 75 * to this distribution, this method returns P(x0 ≤ X ≤ x1). 76 * 77 * @param x0 the (inclusive) lower bound 78 * @param x1 the (inclusive) upper bound 79 * @return the probability that a random variable with this distribution 80 * will take a value between <code>x0</code> and <code>x1</code>, 81 * including the endpoints. 82 * @throws MathException if the cumulative probability can not be 83 * computed due to convergence or other numerical errors. 84 * @throws IllegalArgumentException if <code>x0 > x1</code> 85 */ 86 @Override 87 public double cumulativeProbability(double x0, double x1) 88 throws MathException { 89 if (x0 > x1) { 90 throw MathRuntimeException.createIllegalArgumentException( 91 LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT, x0, x1); 92 } 93 if (FastMath.floor(x0) < x0) { 94 return cumulativeProbability(((int) FastMath.floor(x0)) + 1, 95 (int) FastMath.floor(x1)); // don't want to count mass below x0 96 } else { // x0 is mathematical integer, so use as is 97 return cumulativeProbability((int) FastMath.floor(x0), 98 (int) FastMath.floor(x1)); 99 } 100 } 101 102 /** 103 * For a random variable X whose values are distributed according 104 * to this distribution, this method returns P(X ≤ x). In other words, 105 * this method represents the probability distribution function, or PDF, 106 * for this distribution. 107 * 108 * @param x the value at which the PDF is evaluated. 109 * @return PDF for this distribution. 110 * @throws MathException if the cumulative probability can not be 111 * computed due to convergence or other numerical errors. 112 */ 113 public abstract double cumulativeProbability(int x) throws MathException; 114 115 /** 116 * For a random variable X whose values are distributed according 117 * to this distribution, this method returns P(X = x). In other words, this 118 * method represents the probability mass function, or PMF, for the distribution. 119 * <p> 120 * If <code>x</code> does not represent an integer value, 0 is returned. 121 * 122 * @param x the value at which the probability density function is evaluated 123 * @return the value of the probability density function at x 124 */ 125 public double probability(double x) { 126 double fl = FastMath.floor(x); 127 if (fl == x) { 128 return this.probability((int) x); 129 } else { 130 return 0; 131 } 132 } 133 134 /** 135 * For a random variable X whose values are distributed according 136 * to this distribution, this method returns P(x0 ≤ X ≤ x1). 137 * 138 * @param x0 the inclusive, lower bound 139 * @param x1 the inclusive, upper bound 140 * @return the cumulative probability. 141 * @throws MathException if the cumulative probability can not be 142 * computed due to convergence or other numerical errors. 143 * @throws IllegalArgumentException if x0 > x1 144 */ 145 public double cumulativeProbability(int x0, int x1) throws MathException { 146 if (x0 > x1) { 147 throw MathRuntimeException.createIllegalArgumentException( 148 LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT, x0, x1); 149 } 150 return cumulativeProbability(x1) - cumulativeProbability(x0 - 1); 151 } 152 153 /** 154 * For a random variable X whose values are distributed according 155 * to this distribution, this method returns the largest x, such 156 * that P(X ≤ x) ≤ <code>p</code>. 157 * 158 * @param p the desired probability 159 * @return the largest x such that P(X ≤ x) <= p 160 * @throws MathException if the inverse cumulative probability can not be 161 * computed due to convergence or other numerical errors. 162 * @throws IllegalArgumentException if p < 0 or p > 1 163 */ 164 public int inverseCumulativeProbability(final double p) throws MathException{ 165 if (p < 0.0 || p > 1.0) { 166 throw MathRuntimeException.createIllegalArgumentException( 167 LocalizedFormats.OUT_OF_RANGE_SIMPLE, p, 0.0, 1.0); 168 } 169 170 // by default, do simple bisection. 171 // subclasses can override if there is a better method. 172 int x0 = getDomainLowerBound(p); 173 int x1 = getDomainUpperBound(p); 174 double pm; 175 while (x0 < x1) { 176 int xm = x0 + (x1 - x0) / 2; 177 pm = checkedCumulativeProbability(xm); 178 if (pm > p) { 179 // update x1 180 if (xm == x1) { 181 // this can happen with integer division 182 // simply decrement x1 183 --x1; 184 } else { 185 // update x1 normally 186 x1 = xm; 187 } 188 } else { 189 // update x0 190 if (xm == x0) { 191 // this can happen with integer division 192 // simply increment x0 193 ++x0; 194 } else { 195 // update x0 normally 196 x0 = xm; 197 } 198 } 199 } 200 201 // insure x0 is the correct critical point 202 pm = checkedCumulativeProbability(x0); 203 while (pm > p) { 204 --x0; 205 pm = checkedCumulativeProbability(x0); 206 } 207 208 return x0; 209 } 210 211 /** 212 * Reseeds the random generator used to generate samples. 213 * 214 * @param seed the new seed 215 * @since 2.2 216 */ 217 public void reseedRandomGenerator(long seed) { 218 randomData.reSeed(seed); 219 } 220 221 /** 222 * Generates a random value sampled from this distribution. The default 223 * implementation uses the 224 * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling"> inversion method.</a> 225 * 226 * @return random value 227 * @since 2.2 228 * @throws MathException if an error occurs generating the random value 229 */ 230 public int sample() throws MathException { 231 return randomData.nextInversionDeviate(this); 232 } 233 234 /** 235 * Generates a random sample from the distribution. The default implementation 236 * generates the sample by calling {@link #sample()} in a loop. 237 * 238 * @param sampleSize number of random values to generate 239 * @since 2.2 240 * @return an array representing the random sample 241 * @throws MathException if an error occurs generating the sample 242 * @throws IllegalArgumentException if sampleSize is not positive 243 */ 244 public int[] sample(int sampleSize) throws MathException { 245 if (sampleSize <= 0) { 246 MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NOT_POSITIVE_SAMPLE_SIZE, sampleSize); 247 } 248 int[] out = new int[sampleSize]; 249 for (int i = 0; i < sampleSize; i++) { 250 out[i] = sample(); 251 } 252 return out; 253 } 254 255 /** 256 * Computes the cumulative probability function and checks for NaN values returned. 257 * Throws MathException if the value is NaN. Rethrows any MathException encountered 258 * evaluating the cumulative probability function. Throws 259 * MathException if the cumulative probability function returns NaN. 260 * 261 * @param argument input value 262 * @return cumulative probability 263 * @throws MathException if the cumulative probability is NaN 264 */ 265 private double checkedCumulativeProbability(int argument) throws MathException { 266 double result = Double.NaN; 267 result = cumulativeProbability(argument); 268 if (Double.isNaN(result)) { 269 throw new MathException(LocalizedFormats.DISCRETE_CUMULATIVE_PROBABILITY_RETURNED_NAN, argument); 270 } 271 return result; 272 } 273 274 /** 275 * Access the domain value lower bound, based on <code>p</code>, used to 276 * bracket a PDF root. This method is used by 277 * {@link #inverseCumulativeProbability(double)} to find critical values. 278 * 279 * @param p the desired probability for the critical value 280 * @return domain value lower bound, i.e. 281 * P(X < <i>lower bound</i>) < <code>p</code> 282 */ 283 protected abstract int getDomainLowerBound(double p); 284 285 /** 286 * Access the domain value upper bound, based on <code>p</code>, used to 287 * bracket a PDF root. This method is used by 288 * {@link #inverseCumulativeProbability(double)} to find critical values. 289 * 290 * @param p the desired probability for the critical value 291 * @return domain value upper bound, i.e. 292 * P(X < <i>upper bound</i>) > <code>p</code> 293 */ 294 protected abstract int getDomainUpperBound(double p); 295 296 /** 297 * Use this method to get information about whether the lower bound 298 * of the support is inclusive or not. For discrete support, 299 * only true here is meaningful. 300 * 301 * @return true (always but at Integer.MIN_VALUE because of the nature of discrete support) 302 * @since 2.2 303 */ 304 public boolean isSupportLowerBoundInclusive() { 305 return true; 306 } 307 308 /** 309 * Use this method to get information about whether the upper bound 310 * of the support is inclusive or not. For discrete support, 311 * only true here is meaningful. 312 * 313 * @return true (always but at Integer.MAX_VALUE because of the nature of discrete support) 314 * @since 2.2 315 */ 316 public boolean isSupportUpperBoundInclusive() { 317 return true; 318 } 319 } 320
