rvms.java

Simulation Modeling,Statistical Analysis of Simulation Models,Discrete Event Simulation

/* ------------------------------------------------------------------------- 
 * This is a Java library that can be used to evaluate the probability 
 * density functions (pdf's), cumulative distribution functions (cdf's), and 
 * inverse distribution functions (idf's) for a variety of discrete and 
 * continuous random variables.
 *
 * The following notational conventions are used
 *                 x : possible value of the random variable
 *                 u : real variable (probability) between 0.0 and 1.0 
 *  a, b, n, p, m, s : distribution-specific parameters
 *
 * There are pdf's, cdf's and idf's for 6 discrete random variables
 *
 *      Random Variable    Range (x)  Mean         Variance
 *
 *      bernoulli(p)       0..1       p            p*(1-p)
 *      binomial(n, p)     0..n       n*p          n*p*(1-p)
 *      equilikely(a, b)   a..b       (a+b)/2      ((b-a+1)*(b-a+1)-1)/12 
 *      geometric(p)       0...       p/(1-p)      p/((1-p)*(1-p))
 *      pascal(n, p)       0...       n*p/(1-p)    n*p/((1-p)*(1-p))
 *      poisson(m)         0...       m            m
 *
 * and for 7 continuous random variables
 *
 *      uniform(a, b)      a < x < b  (a+b)/2      (b-a)*(b-a)/12
 *      exponential(m)     x > 0      m            m*m
 *      erlang(n, b)       x > 0      n*b          n*b*b
 *      normal(m, s)       all x      m            s*s
 *      logNormal(a, b)    x > 0         see below
 *      chiSquare(n)       x > 0      n            2*n
 *      student(n)         all x      0  (n > 1)   n/(n-2)   (n > 2)
 *
 * For the Lognormal(a, b), the mean and variance are
 *
 *                        mean = Exp(a + 0.5*b*b)
 *                    variance = (Exp(b*b) - 1)*Exp(2*a + b*b)
 *
 * Name            : Rvms.java (Random Variable ModelS)
 * Authors         : Steve Park & Dave Geyer
 * Translated by   : Richard Dutton & Jun Wang
 * Language        : Java
 * Latest Revision : 7-1-04
 * ------------------------------------------------------------------------- 
 */

import java.math.*;

public class Rvms{
    final double TINY = 1.0e-10;
    final double SQRT2PI = 2.506628274631;	/* sqrt(2 * pi) */
    
    public Rvms(){}
    

    public double pdfBernoulli(double p, long x)
/* =======================================
 * NOTE: use 0.0 < p < 1.0 and 0 <= x <= 1
 * =======================================
 */
    {
	return ((x == 0) ? 1.0 - p : p);
    }

    public double cdfBernoulli(double p, long x)
/* =======================================
 * NOTE: use 0.0 < p < 1.0 and 0 <= x <= 1 
 * =======================================
 */
    {
	return ((x == 0) ? 1.0 - p : 1.0);
    }

    public long idfBernoulli(double p, double u)
/* =========================================
 * NOTE: use 0.0 < p < 1.0 and 0.0 < u < 1.0 
 * =========================================
 */
    {
	return ((u < 1.0 - p) ? 0 : 1);
    }

    public double pdfEquilikely(long a, long b, long x)
/* ============================================ 
 * NOTE: use a <= x <= b 
 * ============================================
 */
    {
	return (1.0 / (b - a + 1.0));
    }

    public double cdfEquilikely(long a, long b, long x)
/* ============================================
 * NOTE: use a <= x <= b 
 * ============================================
 */
    {
	return ((x - a + 1.0) / (b - a + 1.0));
    }

    public long idfEquilikely(long a, long b, double u)
/* ============================================ 
 * NOTE: use a <= b and 0.0 < u < 1.0 
 * ============================================
 */
    {
	return (a + (long) (u * (b - a + 1)));
    }

    public double pdfBinomial(long n, double p, long x)
/* ============================================ 
 * NOTE: use 0 <= x <= n and 0.0 < p < 1.0 
 * ============================================
 */
    {
	double s, t;
	
	s = logChoose(n, x);
	t = x * Math.log(p) + (n - x) * Math.log(1.0 - p);
	return (Math.exp(s + t));
    }

    public double cdfBinomial(long n, double p, long x)
/* ============================================ 
 * NOTE: use 0 <= x <= n and 0.0 < p < 1.0 
 * ============================================
 */
    {
	if (x < n)
	    return (1.0 - inBeta(x + 1, n - x, p));
	else
	    return (1.0);
    }

    public long idfBinomial(long n, double p, double u)
/* ================================================= 
 * NOTE: use 0 <= n, 0.0 < p < 1.0 and 0.0 < u < 1.0 
 * =================================================
 */
    {
	long x = (long) (n * p);             /* start searching at the mean */
	
	if (cdfBinomial(n, p, x) <= u)
	    while (cdfBinomial(n, p, x) <= u)
		x++;
	else if (cdfBinomial(n, p, 0) <= u)
	    while (cdfBinomial(n, p, x - 1) > u)
		x--;
	else
	    x = 0;
	return (x);
    }

    public double pdfGeometric(double p, long x)
/* ===================================== 
 * NOTE: use 0.0 < p < 1.0 and x >= 0 
 * =====================================
 */
    {
	return ((1.0 - p) * Math.exp(x * Math.log(p)));
    }

    public double cdfGeometric(double p, long x)
/* ===================================== 
 * NOTE: use 0.0 < p < 1.0 and x >= 0 
 * =====================================
 */
    {
	return (1.0 - Math.exp((x + 1) * Math.log(p)));
    }

    public long idfGeometric(double p, double u)
/* ========================================= 
 * NOTE: use 0.0 < p < 1.0 and 0.0 < u < 1.0 
 * =========================================
 */
    {
	return ((long) (Math.log(1.0 - u) / Math.log(p)));
    }

    public double pdfPascal(long n, double p, long x)
/* =========================================== 
 * NOTE: use n >= 1, 0.0 < p < 1.0, and x >= 0 
 * ===========================================
 */
    {
	double  s, t;
	
	s = logChoose(n + x - 1, x);
	t = x * Math.log(p) + n * Math.log(1.0 - p);
	return (Math.exp(s + t));
    }

    public double cdfPascal(long n, double p, long x)
/* =========================================== 
 * NOTE: use n >= 1, 0.0 < p < 1.0, and x >= 0 
 * ===========================================
 */
    {
	return (1.0 - inBeta(x + 1, n, p));
    }

    public long idfPascal(long n, double p, double u)
/* ================================================== 
 * NOTE: use n >= 1, 0.0 < p < 1.0, and 0.0 < u < 1.0 
 * ==================================================
 */
    {
	long x = (long) (n * p / (1.0 - p));    /* start searching at the mean */
	
	if (cdfPascal(n, p, x) <= u)
	    while (cdfPascal(n, p, x) <= u)
		x++;
	else if (cdfPascal(n, p, 0) <= u)
	    while (cdfPascal(n, p, x - 1) > u)
		x--;
	else
	    x = 0;
	return (x);
    }

    public double pdfPoisson(double m, long x)
/* ===================================
 * NOTE: use m > 0 and x >= 0 
 * ===================================
 */
    {
	double t;
	
	t = - m + x * Math.log(m) - logFactorial(x);
	return (Math.exp(t));
    }

    public double cdfPoisson(double m, long x)
/* =================================== 
 * NOTE: use m > 0 and x >= 0 
 * ===================================
 */
    {
	return (1.0 - inGamma(x + 1, m));
    }

    public long idfPoisson(double m, double u)
/* =================================== 
 * NOTE: use m > 0 and 0.0 < u < 1.0 
 * ===================================
 */
    {
	long x = (long) m;                    /* start searching at the mean */
	
	if (cdfPoisson(m, x) <= u)
	    while (cdfPoisson(m, x) <= u)
		x++;
	else if (cdfPoisson(m, 0) <= u)
	    while (cdfPoisson(m, x - 1) > u)
		x--;
	else
	    x = 0;
	return (x);
    }

    public double pdfUniform(double a, double b, double x)
/* =============================================== 
 * NOTE: use a < x < b 
 * ===============================================
 */
    {
	return (1.0 / (b - a));
    }

    public double cdfUniform(double a, double b, double x)
/* =============================================== 
 * NOTE: use a < x < b 
 * ===============================================
 */
    {
	return ((x - a) / (b - a));
    }

    public double idfUniform(double a, double b, double u)
/* =============================================== 
 * NOTE: use a < b and 0.0 < u < 1.0 
 * ===============================================
 */
    {
	return (a + (b - a) * u);
    }

    public double pdfExponential(double m, double x)
/* ========================================= 
 * NOTE: use m > 0 and x > 0 
 * =========================================
 */
    {
	return ((1.0 / m) * Math.exp(- x / m));
    }

    public double cdfExponential(double m, double x)
/* ========================================= 
 * NOTE: use m > 0 and x > 0 
 * =========================================
 */
    {
	return (1.0 - Math.exp(- x / m));
    }

    public double idfExponential(double m, double u)
/* ========================================= 
 * NOTE: use m > 0 and 0.0 < u < 1.0 
 * =========================================
 */
    {
	return (- m * Math.log(1.0 - u));
    }

    public double pdfErlang(long n, double b, double x)
	/* ============================================ 
	 * NOTE: use n >= 1, b > 0, and x > 0 
	 * ============================================
	 */
    {
	double t;

	t = (n - 1) * Math.log(x / b) - (x / b) - Math.log(b) - logGamma(n);
	return (Math.exp(t));
    }

    public double cdfErlang(long n, double b, double x)
	/* ============================================ 
	 * NOTE: use n >= 1, b > 0, and x > 0 
	 * ============================================
	 */
    {
	return (inGamma(n, x / b));
    }

    public double idfErlang(long n, double b, double u)
	/* ============================================ 
	 * NOTE: use n >= 1, b > 0 and 0.0 < u < 1.0 
	 * ============================================
	 */
    {
	double t, x = n * b;                   /* initialize to the mean, then */

	do {                                   /* use Newton-Raphson iteration */
	    t = x;
	    x = t + (u - cdfErlang(n, b, t)) / pdfErlang(n, b, t);
	    if (x <= 0.0)
		x = 0.5 * t;
	} while (Math.abs(x - t) >= TINY);
	return (x);
    }

    public double pdfStandard(double x)
	/* =================================== 
	 * NOTE: x can be any value 
	 * ===================================
	 */
    {
	return (Math.exp(- 0.5 * x * x) / SQRT2PI);
    }

    public double cdfStandard(double x)
	/* =================================== 
	 * NOTE: x can be any value 
	 * ===================================
	 */
    { 
	double t;

	t = inGamma(0.5, 0.5 * x * x);
	if (x < 0.0)
	    return (0.5 * (1.0 - t));
	else