/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
package umontreal.iro.lecuyer.probdist;
import umontreal.iro.lecuyer.util.*;
import umontreal.iro.lecuyer.functions.MathFunction;
import optimization.*;
/**
* Extends the class {@link DiscreteDistributionInt} for
* the negative binomial distribution with real
* parameters γ and p, where
* γ > 0 and
* 0 <= p <= 1.
* Its mass function is
*
*
* where Γ is the gamma function.
*
*
* If γ is an integer, p(x) can be interpreted as the probability
* of having x failures before the γ-th success in a sequence of
* independent Bernoulli trials with probability of success p. This special
* case is implemented as the Pascal distribution (see {@link PascalDist}).
*
*/
public class NegativeBinomialDist extends DiscreteDistributionInt {
protected double gamma;
protected double p;
private static final double EPSI = 1.0E-10;
private static class Func1 implements MathFunction {
protected int n;
protected int[] x;
protected double p;
public Func1 (double p, int[] x, int n) {
this.p = p;
this.n = n;
this.x = x;
}
public double evaluate (double gam) {
if (gam <= 0 ) return 1.0e100;
double sum = 0.0;
for (int j = 0; j < n; j++)
sum += Num.digamma (gam + x[j]);
return sum/n + Math.log (p) - Num.digamma (gam);
}
}
private static class Function implements MathFunction {
protected int m;
protected int max;
protected double mean;
protected int[] Fj;
public Function (int m, int max, double mean, int[] Fj) {
this.m = m;
this.max = max;
this.mean = mean;
this.Fj = new int[Fj.length];
System.arraycopy(Fj, 0, this.Fj, 0, Fj.length);
}
public double evaluate (double s) {
if (s <= 0 ) return 1.0e100;
double sum = 0.0;
double p = s / (s + mean);
for (int j = 0; j < max; j++)
sum += Fj[j] / (s + (double) j);
return sum + m * Math.log (p);
}
}
/*
// Class Optim seems to be useless
private static class Optim implements Lmder_fcn
{
private double mean;
private int N;
private int max;
private int [] Fj;
public Optim (int N, int max, double mean, int[] Fj)
{
this.N = N;
this.max = max;
this.mean = mean;
this.Fj = new int[max];
System.arraycopy (Fj, 0, this.Fj, 0, max);
}
public void fcn (int m, int n, double[] x, double[] fvec, double[][] fjac,
int iflag[])
{
if (x[1] <= 0.0 || x[2] <= 0.0 || x[2] >= 1.0) {
final double BIG = 1.0e100;
fvec[1] = BIG;
fvec[2] = BIG;
fjac[1][1] = BIG;
fjac[1][2] = 0.0;
fjac[2][1] = 0.0;
fjac[2][2] = BIG;
return;
}
double trig;
if (iflag[1] == 1)
{
double sum = 0.0;
for (int j = 0; j < max; j++)
sum += Fj[j] / (x[1] + j);
fvec[1] = x[2] * mean - x[1] * (1.0 - x[2]);
fvec[2] = N * Math.log (x[2]) + sum;
} else if (iflag[1] == 2) {
fjac[1][1] = x[2] - 1.0;
fjac[1][2] = mean + x[1];
double sum = 0.0;
for (int j = 0; j < max; j++)
sum += Fj[j] / ((x[1] + j)*(x[1] + j));
fjac[2][1] = -sum;
fjac[2][2] = N / x[2];
}
}
}
*/
public static double MAXN = 100000;
protected NegativeBinomialDist() {}
/**
* Creates an object that contains the probability
* terms and the distribution function for
* the negative binomial distribution with parameters γ and p.
*
*/
public NegativeBinomialDist (double gamma, double p) {
setParams (gamma, p);
}
public double prob (int x) {
if (x < 0)
return 0.0;
if (p == 0.0)
return 0.0;
if (p == 1.0) {
if (x > 0)
return 0.0;
else
return 1.0;
}
if (pdf == null)
return prob (gamma, p, x);
if (x > xmax || x < xmin)
return prob (gamma, p, x);
return pdf[x - xmin];
}
public double cdf (int x) {
if (x < 0)
return 0.0;
if (p >= 1.0) // In fact, p == 1
return 1.0;
if (p <= 0.0) // In fact, p == 0
return 0.0;
if (cdf != null) {
if (x >= xmax)
return 1.0;
if (x < xmin)
return cdf (gamma, p, x);
if (x <= xmed)
return cdf[x - xmin];
else
// We keep the complementary distribution in the upper part of cdf
return 1.0 - cdf[x + 1 - xmin];
}
else
return cdf (gamma, p, x);
}
public double barF (int x) {
if (x < 1)
return 1.0;
if (p >= 1.0) // In fact, p == 1
return 0.0;
if (p <= 0.0) // In fact, p == 0
return 1.0;
if (cdf == null)
//return BinomialDist.cdf (x - 1 + n, p, n - 1);
return BetaDist.barF (gamma, x, 15, p);
if (x > xmax)
//return BinomialDist.cdf (x - 1 + n, p, n - 1);
return BetaDist.barF (gamma, x, 15, p);
if (x <= xmin)
return 1.0;
if (x > xmed)
// We keep the complementary distribution in the upper part of cdf
return cdf[x - xmin];
else
return 1.0 - cdf[x - 1 - xmin];
}
public int inverseFInt (double u) {
if (cdf == null)
return inverseF (gamma, p, u);
else
return super.inverseFInt (u);
}
public double getMean() {
return NegativeBinomialDist.getMean (gamma, p);
}
public double getVariance() {
return NegativeBinomialDist.getVariance (gamma, p);
}
public double getStandardDeviation() {
return NegativeBinomialDist.getStandardDeviation (gamma, p);
}
/**
* Computes the probability p(x).
*
*/
public static double prob (double gamma, double p, int x) {
final int SLIM = 15; // To avoid overflow
final double MAXEXP = (Num.DBL_MAX_EXP - 1)*Num.LN2;// To avoid overflow
final double MINEXP = (Num.DBL_MIN_EXP - 1)*Num.LN2;// To avoid underflow
double y;
if (p < 0.0 || p > 1.0)
throw new IllegalArgumentException ("p not in [0, 1]");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0.0");
if (x < 0)
return 0.0;
if (p >= 1.0) { // In fact, p == 1
if (x == 0)
return 1.0;
else
return 0.0;
}
if (p <= 0.0) // In fact, p == 0
return 0.0;
y = Num.lnGamma (gamma + x) - (Num.lnFactorial (x) + Num.lnGamma (gamma))
+ gamma * Math.log (p) + x * Math.log1p (-p) ;
if (y >= MAXEXP)
throw new IllegalArgumentException ("term overflow");
if (y <= MINEXP)
return 0.0;
else
return Math.exp (y);
}
/**
* Computes the distribution function.
*
*/
public static double cdf (double gamma, double p, int x) {
final double EPSILON = DiscreteDistributionInt.EPSILON;
final int LIM1 = 100000;
double sum, term, termmode;
int i, mode;
if (p < 0.0 || p > 1.0)
throw new IllegalArgumentException ("p not in [0, 1]");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0.0");
if (x < 0)
return 0.0;
if (p >= 1.0) // In fact, p == 1
return 1.0;
if (p <= 0.0) // In fact, p == 0
return 0.0;
// Compute the maximum term
mode = 1 + (int) Math.floor ((gamma*(1.0 - p) - 1.0)/p);
if (mode < 0)
mode = 0;
else if (mode > x)
mode = x;
if (mode <= LIM1) {
sum = term = termmode = prob (gamma, p, mode);
for (i = mode; i > 0; i--) {
term *= i/((1.0 - p)*(gamma + i - 1.0));
if (term < EPSILON)
break;
sum += term;
}
term = termmode;
for (i = mode; i < x; i++) {
term *= (1.0 - p)*(gamma + i)/(i + 1);
if (term < EPSILON)
break;
sum += term;
}
if (sum <= 1.0)
return sum;
else
return 1.0;
}
else
//return 1.0 - BinomialDist.cdf (x + n, p, n - 1);
return BetaDist.cdf (gamma, x + 1.0, 15, p);
}
/**
* Computes the inverse function without precomputing tables.
* This method computes the CDF at the mode (maximum term) and performs
* a linear search from that point.
*
*/
public static int inverseF (double gamma, double p, double u) {
// final int LIM1 = Integer.MAX_VALUE; //100000;
if (u < 0.0 || u >= 1.0)
throw new IllegalArgumentException ("u is not in [0,1]");
if (p < 0.0 || p > 1.0)
throw new IllegalArgumentException ("p not in [0, 1]");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0");
if (p >= 1.0) // In fact, p == 1
return 0;
if (p <= 0.0) // In fact, p == 0
return 0;
if (u <= 0.0)
return 0;
if (u >= 1.0)
return Integer.MAX_VALUE;
double sum, term, termmode;
int i, mode, x = 0;
// Compute the maximum term
mode = 1 + (int) Math.floor ((gamma*(1.0 - p) - 1.0)/p);
if (mode < 0)
mode = 0;
// if (mode <= LIM1) { // This test does not seem necessary
// We compute the CDF at the mode.
sum = term = termmode = prob (gamma, p, mode);
for (i = mode; i > 0; i--) {
term *= i/((1.0 - p)*(gamma + i - 1.0));
if (term < EPSILON)
break;
sum += term;
}
term = termmode;
x = mode;
if (sum < u) {
// The CDF at the mode is less than u, so we add term to get >= u.
while (sum < u) {
term *= (1.0 - p)*(gamma + x)/(x + 1);
if (term < EPSILON)
break;
sum += term;
x++;
}
}
else {
// The computed CDF is too big so we substract from it.
sum -= term;
while (sum >= u) {
term *= x/((1.0 - p)*(gamma + x - 1.0));
x--;
if (term < EPSILON)
break;
sum -= term;
}
}
return x;
}
/**
* Estimates the parameter p of the negative binomial distribution
* using the maximum likelihood method, from the n observations
* x[i],
* i = 0, 1,…, n - 1. The parameter
* γ = gamma is assumed known.
* The estimate hat(p) is returned in element 0
* of the returned array.
* The maximum likelihood estimator hat(p) satisfies the equation
*
* hat(p) = γ/(γ + bar(x)n),
* where bar(x)n is the average of
* x[0],…, x[n - 1].
*
* @param x the list of observations used to evaluate parameters
*
* @param n the number of observations used to evaluate parameters
*
* @param gamma the first parameter of the negative binomial
*
* @return returns the parameters [hat(p)]
*
*/
public static double[] getMLE (int[] x, int n, double gamma) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
double mean = 0.0;
for (int i = 0; i < n; i++) {
mean += x[i];
}
mean /= (double) n;
double[] param = new double[1];
param[0] = gamma / (gamma + mean);
return param;
}
/**
* Creates a new instance of a negative binomial distribution with parameters
* γ = gamma given and hat(p) estimated using the maximum
* likelihood method, from the n observations x[i],
*
* i = 0, 1,…, n - 1.
*
* @param x the list of observations to use to evaluate parameters
*
* @param n the number of observations to use to evaluate parameters
*
* @param gamma the first parameter of the negative binomial
*
*
*/
public static NegativeBinomialDist getInstanceFromMLE (int[] x, int n,
double gamma) {
double parameters[] = getMLE (x, n, gamma);
return new NegativeBinomialDist (gamma, parameters[0]);
}
/**
* Estimates the parameter γ of the negative binomial distribution
* using the maximum likelihood method, from the n observations
* x[i],
* i = 0, 1,…, n - 1. The parameter p is assumed known.
* The estimate
* hat(γ) is returned in element 0 of the returned array.
*
* @param x the list of observations used to evaluate parameters
*
* @param n the number of observations used to evaluate parameters
*
* @param p the second parameter of the negative binomial
*
* @return returns the parameters [
* hat(γ)]
*
*/
public static double[] getMLE1 (int[] x, int n, double p) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
double mean = 0.0;
for (int i = 0; i < n; i++)
mean += x[i];
mean /= n;
double gam0 = mean*p/(1.0 - p);
double[] param = new double[1];
Func1 f = new Func1 (p, x, n);
param[0] = RootFinder.brentDekker (gam0/10.0, 10.0*gam0, f, 1e-5);
return param;
}
/**
* Creates a new instance of a negative binomial distribution with parameters
* p given and
* hat(γ) estimated using the maximum
* likelihood method, from the n observations x[i],
*
* i = 0, 1,…, n - 1.
*
* @param x the list of observations to use to evaluate parameters
*
* @param n the number of observations to use to evaluate parameters
*
* @param p the second parameter of the negative binomial
*
*
*/
public static NegativeBinomialDist getInstanceFromMLE1 (int[] x, int n,
double p) {
double param[] = getMLE1 (x, n, p);
return new NegativeBinomialDist (param[0], p);
}
/**
* Estimates the parameter
* (γ, p) of the negative binomial distribution
* using the maximum likelihood method, from the n observations
* x[i],
* i = 0, 1,…, n - 1. The estimates are returned in a two-element
* array, in regular order: [γ, p].
*
* @param x the list of observations used to evaluate parameters
*
* @param n the number of observations used to evaluate parameters
*
* @return returns the parameters [
* hat(γ), hat(p)]
*
*/
public static double[] getMLE (int[] x, int n) {
double estimGamma;
if (n <= 0)
throw new IllegalArgumentException ("n<= 0");
double sum = 0.0;
double max = Integer.MIN_VALUE;
for (int i = 0; i < n; i++)
{
sum += x[i];
if (x[i] > max)
max = x[i];
}
double mean = (double) sum / (double) n;
double var = 0.0;
for (int i = 0; i < n; i++)
var += (x[i] - mean) * (x[i] - mean);
var /= (double) n;
if (mean >= var)
throw new UnsupportedOperationException("mean >= variance");
estimGamma = (mean * mean) / ( var - mean );
int[] Fj = new int[(int) max];
for (int j = 0; j < max; j++)
{
int prop = 0;
for (int i = 0; i < n; i++)
if (x[i] > j)
prop++;
Fj[j] = prop;
}
double[] param = new double[3];
Function f = new Function (n, (int) max, mean, Fj);
param[1] = RootFinder.brentDekker (estimGamma/10, estimGamma*10, f, 1e-5);
param[2] = param[1] / (param[1] + mean);
/* Seems to be useless
Optim system = new Optim (n, (int) max, mean, Fj);
double[] fvec = new double [3];
double[][] fjac = new double[3][3];
int[] iflag = new int[2];
int[] info = new int[2];
int[] ipvt = new int[3];
Minpack_f77.lmder1_f77 (system, 2, 2, param, fvec, fjac, 1e-5, info, ipvt);
*/
double parameters[] = new double[2];
parameters[0] = param[1];
parameters[1] = param[2];
return parameters;
}
/**
* Creates a new instance of a negative binomial distribution with
* parameters γ and p estimated using the maximum likelihood method
* based on the n observations x[i],
* i = 0, 1,…, n - 1.
*
* @param x the list of observations to use to evaluate parameters
*
* @param n the number of observations to use to evaluate parameters
*
*
*/
public static NegativeBinomialDist getInstanceFromMLE (int[] x, int n) {
double parameters[] = getMLE (x, n);
return new NegativeBinomialDist (parameters[0], parameters[1]);
}
/**
* Computes and returns the mean
* E[X] = γ(1 - p)/p
* of the negative binomial distribution with parameters γ and p.
*
* @return the mean of the negative binomial distribution
*
* E[X] = γ(1 - p)/p
*
*/
public static double getMean (double gamma, double p) {
if (p < 0.0 || p > 1.0)
throw new IllegalArgumentException ("p not in [0, 1]");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0");
return (gamma * (1 - p) / p);
}
/**
* Computes and returns the variance
* Var[X] = γ(1 - p)/p2
* of the negative binomial distribution with parameters γ and p.
*
* @return the variance of the negative binomial distribution
*
* Var[X] = γ(1 - p)/p2
*
*/
public static double getVariance (double gamma, double p) {
if (p < 0.0 || p > 1.0)
throw new IllegalArgumentException ("p not in [0, 1]");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0");
return (gamma * (1 - p) / (p * p));
}
/**
* Computes and returns the standard deviation of the negative
* binomial distribution with parameters γ and p.
*
* @return the standard deviation of the negative binomial distribution
*
*/
public static double getStandardDeviation (double gamma, double p) {
return Math.sqrt (NegativeBinomialDist.getVariance (gamma, p));
}
/**
* Returns the parameter γ of this object.
*
*/
public double getGamma() {
return gamma;
}
/**
* Returns the parameter p of this object.
*
*/
public double getP() {
return p;
}
/**
* Sets the parameter γ and p of this object.
*
*/
public void setParams (double gamma, double p) {
/* *
* Compute all probability terms of the negative binomial distribution;
* start at the mode, and calculate probabilities on each side until they
* become smaller than EPSILON. Set all others to 0.
*/
supportA = 0;
int i, mode, Nmax;
int imin, imax;
double sum;
double[] P; // Negative Binomial mass probabilities
double[] F; // Negative Binomial cumulative
if (p < 0.0 || p > 1.0)
throw new IllegalArgumentException ("p not in [0, 1]");
if (gamma <= 0.0)
throw new IllegalArgumentException ("gamma <= 0");
this.gamma = gamma;
this.p = p;
// Compute the mode (at the maximum term)
mode = 1 + (int) Math.floor((gamma*(1.0 - p) - 1.0)/p);
/* *
For mode > MAXN, we shall not use pre-computed arrays.
mode < 0 should be impossible, unless overflow of long occur, in
which case mode will be = LONG_MIN.
*/
if (mode < 0.0 || mode > MAXN) {
pdf = null;
cdf = null;
return;
}
/* *
In theory, the negative binomial distribution has an infinite range.
But for i > Nmax, probabilities should be extremely small.
Nmax = Mean + 16 * Standard deviation.
*/
Nmax = (int)(gamma*(1.0 - p)/p + 16*Math.sqrt (gamma*(1.0 - p)/(p*p)));
if (Nmax < 32)
Nmax = 32;
P = new double[1 + Nmax];
F = new double[1 + Nmax];
double epsilon = EPSILON/prob (gamma, p, mode);
// We shall normalize by explicitly summing all terms >= epsilon
sum = P[mode] = 1.0;
// Start from the maximum and compute terms > epsilon on each side.
i = mode;
while (i > 0 && P[i] >= epsilon) {
P[i - 1] = P[i]*i/((1.0 - p)*(gamma + i - 1));
i--;
sum += P[i];
}
imin = i;
i = mode;
while (P[i] >= epsilon) {
P[i + 1] = P[i]*(1.0 - p)*(gamma + i)/(i + 1);
i++;
sum += P[i];
if (i == Nmax - 1) {
Nmax *= 2;
double[] nT = new double[1 + Nmax];
System.arraycopy (P, 0, nT, 0, P.length);
P = nT;
nT = new double[1 + Nmax];
System.arraycopy (F, 0, nT, 0, F.length);
F = nT;
}
}
imax = i;
// Renormalize the sum of probabilities to 1
for (i = imin; i <= imax; i++)
P[i] /= sum;
// Compute the cumulative probabilities for F and keep them in the
// lower part of CDF.
F[imin] = P[imin];
i = imin;
while (i < imax && F[i] < 0.5) {
i++;
F[i] = F[i - 1] + P[i];
}
// This is the boundary between F (i <= xmed) and 1 - F (i > xmed) in
// the array CDF
xmed = i;
// Compute the cumulative probabilities of the complementary
// distribution 1 - F and keep them in the upper part of the array
F[imax] = P[imax];
i = imax - 1;
do {
F[i] = P[i] + F[i + 1];
i--;
} while (i > xmed);
xmin = imin;
xmax = imax;
pdf = new double[imax + 1 - imin];
cdf = new double[imax + 1 - imin];
System.arraycopy (P, imin, pdf, 0, imax + 1 - imin);
System.arraycopy (F, imin, cdf, 0, imax + 1 - imin);
}
/**
* Return a table containing the parameters of the current distribution.
* This table is put in regular order: [γ, p].
*
*
*/
public double[] getParams () {
double[] retour = {gamma, p};
return retour;
}
public String toString () {
return getClass().getSimpleName() + " : gamma = " + gamma + ", p = " + p;
}
}
Ads help maintain this website.