package umontreal.iro.lecuyer.probdist;
/**
* Classes implementing discrete distributions over the integers should
* inherit from this class.
* It specifies the signatures of methods for computing the mass function
* (or probability)
* p(x) = P[X = x], distribution function F(x),
* complementary distribution function bar(F)(x),
* and inverse distribution function F-1(u), for a random variable X
* with a discrete distribution over the integers.
*
*
* WARNING: the complementary distribution function is defined as
*
* bar(F)(j) = P[X >= j] (for integers j, so that for discrete distributions
* in SSJ,
* F(j) + bar(F)(j)≠1 since both include the term P[X = j].
*
*
* The implementing classes provide both static and non-static methods
* to compute the above functions.
* The non-static methods require the creation of an object of
* class {@link umontreal.iro.lecuyer.probdist.DiscreteDistributionInt DiscreteDistributionInt};
* all the non-negligible terms of the mass and distribution functions will be
* precomputed by the constructor and kept in arrays. Subsequent accesses
* will be very fast.
* The static methods do not require the construction of an object.
* These static methods are not specified in this abstract class because
* the number and types of their parameters depend on the distribution.
* When methods have to be called several times
* with the same parameters for the distributions,
* it is usually more efficient to create an object and use its non-static
* methods instead of the static ones.
* This trades memory for speed.
*
*/
public abstract class DiscreteDistributionInt implements Distribution {
/**
* Environment variable that determines what probability terms can
* be considered as negligible when building precomputed tables for
* distribution and mass functions. Probabilities smaller than EPSILON
* are not stored in the
* {@link umontreal.iro.lecuyer.probdist.DiscreteDistribution DiscreteDistribution} objects
* (such as those of class {@link PoissonDist}, etc.), but are computed
* directly each time they are needed (which should be very seldom).
* The default value is set to 10-16.
*
*/
public static double EPSILON = 1.0e-16;
/*
For better precision in the tails, we keep the cumulative probabilities
(F) in cdf[x] for x <= xmed (i.e. cdf[x] is the sum off all the probabi-
lities pdf[i] for i <= x),
and the complementary cumulative probabilities (1 - F) in cdf[x] for
x > xmed (i.e. cdf[x] is the sum off all the probabilities pdf[i]
for i >= x).
*/
protected final static double EPS_EXTRA = 1 / 100.0;
protected double cdf[] = null; // cumulative probabilities
protected double pdf[] = null; // probability terms or mass distribution
protected int xmin = 0; // pdf[x] = 0 for x < xmin
protected int xmax = 0; // pdf[x] = 0 for x > xmax
protected int xmed = 0; // cdf[x] = F (x) for x <= xmed, and
// cdf[x] = bar_F (x) for x > xmed
protected int supportA = Integer.MIN_VALUE;
protected int supportB = Integer.MAX_VALUE;
/**
* Returns p(x), the probability of x,
* which should be a real number in the interval [0, 1].
*
* @param x value at which the mass function must be evaluated
*
* @return the mass function evaluated at x
*
*/
public abstract double prob (int x);
/**
* Returns the distribution function F evaluated at x
* (see).
* Calls the {@link #cdf(int) cdf}(int) method.
*
* @param x value at which the distribution function must be evaluated
*
* @return the distribution function evaluated at x
*
*/
public double cdf (double x) {
return cdf ((int) x);
}
/**
* Returns the distribution function F evaluated at x
* (see).
*
* @param x value at which the distribution function must be evaluated
*
* @return the distribution function evaluated at x
*
*/
public abstract double cdf (int x);
/**
* Returns bar(F)(x), the complementary distribution function.
* Calls the {@link #barF(int) barF}(int) method.
*
* @param x value at which the complementary distribution function
* must be evaluated
*
* @return the complementary distribution function evaluated at x
*
*/
public double barF (double x) {
return barF ((int) x);
}
/**
* Returns bar(F)(x), the complementary
* distribution function. See the WARNING above.
*
* @param x value at which the complementary distribution function
* must be evaluated
*
* @return the complementary distribution function evaluated at x
*
*/
public double barF (int x) {
return 1.0 - cdf (x - 1);
}
/**
* Returns the lower limit xa of the support of the probability
* mass function. The probability is 0 for all x < xa.
*
* @return x lower limit of support
*
*/
public int getXinf() {
return supportA;
}
/**
* Returns the upper limit xb of the support of the probability
* mass function. The probability is 0 for all x > xb.
*
* @return x upper limit of support
*
*/
public int getXsup() {
return supportB;
}
/**
* Returns the inverse distribution function
* F-1(u), where
* 0 <= u <= 1. Calls the inverseFInt method.
*
* @param u value in the interval (0, 1) for which
* the inverse distribution function is evaluated
*
* @return the inverse distribution function evaluated at u
* @exception IllegalArgumentException if u is not in the interval (0, 1)
*
* @exception ArithmeticException if the inverse cannot be computed,
* for example if it would give infinity in a theoritical context
*
*
*/
public double inverseF (double u) {
return inverseFInt (u);
}
/**
* Returns the inverse distribution function
* F-1(u), where
* 0 <= u <= 1.
* The default implementation uses binary search.
*
* @param u value in the interval (0, 1) for which
* the inverse distribution function is evaluated
*
* @return the inverse distribution function evaluated at u
* @exception IllegalArgumentException if u is not in the interval (0, 1)
*
* @exception ArithmeticException if the inverse cannot be computed,
* for example if it would give infinity in a theoritical context
*
*
*/
public int inverseFInt (double u) {
int i, j, k;
if (u < 0.0 || u > 1.0)
throw new IllegalArgumentException ("u is not in [0,1]");
if (u <= 0.0)
return supportA;
if (u >= 1.0)
return supportB;
// Remember: the upper part of cdf contains the complementary distribu-
// tion for xmed < s <= xmax, and the lower part of cdf the
// distribution for xmin <= x <= xmed
if (u <= cdf[xmed - xmin]) {
// In the lower part of cdf
if (u <= cdf[0])
return xmin;
i = 0;
j = xmed - xmin;
while (i < j) {
k = (i + j) / 2;
if (u > cdf[k])
i = k + 1;
else
j = k;
}
}
else {
// In the upper part of cdf
u = 1 - u;
if (u < cdf[xmax - xmin])
return xmax;
i = xmed - xmin + 1;
j = xmax - xmin;
while (i < j) {
k = (i + j) / 2;
if (u < cdf[k])
i = k + 1;
else
j = k;
}
i--;
}
return i + xmin;
}
}
Ads help maintain this website.