/* -------------------------------------------------------------------------  * This is an ANSI C 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.c (Random Variable ModelS) * Author          : Steve Park & Dave Geyer * Language        : ANSI C * Latest Revision : 11-22-97 * -------------------------------------------------------------------------  */#include <math.h>#include "rvms.h"#define TINY    1.0e-10#define SQRT2PI 2.506628274631               /* sqrt(2 * pi) */static double pdfStandard(double x);static double cdfStandard(double x);static double idfStandard(double u);static double LogGamma(double a);static double LogBeta(double a, double b);static double InGamma(double a, double b);static double InBeta(double a, double b, double x);   double pdfBernoulli(double p, long x)/* ======================================= * NOTE: use 0.0 < p < 1.0 and 0 <= x <= 1 * ======================================= */{   return ((x == 0) ? 1.0 - p : p);}   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);}   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);}   double pdfEquilikely(long a, long b, long x)/* ============================================  * NOTE: use a <= x <= b  * ============================================ */{   return (1.0 / (b - a + 1.0));}   double cdfEquilikely(long a, long b, long x)/* ============================================ * NOTE: use a <= x <= b  * ============================================ */{   return ((x - a + 1.0) / (b - a + 1.0));}   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)));}   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 * log(p) + (n - x) * log(1.0 - p);   return (exp(s + t));}   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);}   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);}   double pdfGeometric(double p, long x)/* =====================================  * NOTE: use 0.0 < p < 1.0 and x >= 0  * ===================================== */{   return ((1.0 - p) * exp(x * log(p)));}   double cdfGeometric(double p, long x)/* =====================================  * NOTE: use 0.0 < p < 1.0 and x >= 0  * ===================================== */{   return (1.0 - exp((x + 1) * log(p)));}   long idfGeometric(double p, double u)/* =========================================  * NOTE: use 0.0 < p < 1.0 and 0.0 < u < 1.0  * ========================================= */{   return ((long) (log(1.0 - u) / log(p)));}   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 * log(p) + n * log(1.0 - p);   return (exp(s + t));}   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));}   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);}   double pdfPoisson(double m, long x)/* =================================== * NOTE: use m > 0 and x >= 0  * =================================== */{   double t;   t = - m + x * log(m) - LogFactorial(x);   return (exp(t));}   double cdfPoisson(double m, long x)/* ===================================  * NOTE: use m > 0 and x >= 0  * =================================== */{   return (1.0 - InGamma(x + 1, m));}   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);}   double pdfUniform(double a, double b, double x)/* ===============================================  * NOTE: use a < x < b  * =============================================== */{   return (1.0 / (b - a));}   double cdfUniform(double a, double b, double x)/* ===============================================  * NOTE: use a < x < b  * =============================================== */{   return ((x - a) / (b - a));}   double idfUniform(double a, double b, double u)/* ===============================================  * NOTE: use a < b and 0.0 < u < 1.0  * =============================================== */{   return (a + (b - a) * u);}   double pdfExponential(double m, double x)/* =========================================  * NOTE: use m > 0 and x > 0  * ========================================= */{   return ((1.0 / m) * exp(- x / m));}   double cdfExponential(double m, double x)/* =========================================  * NOTE: use m > 0 and x > 0  * ========================================= */{   return (1.0 - exp(- x / m));}   double idfExponential(double m, double u)/* =========================================  * NOTE: use m > 0 and 0.0 < u < 1.0  * ========================================= */{   return (- m * log(1.0 - u));}   double pdfErlang(long n, double b, double x)/* ============================================  * NOTE: use n >= 1, b > 0, and x > 0  * ============================================ */{   double t;   t = (n - 1) * log(x / b) - (x / b) - log(b) - LogGamma(n);   return (exp(t));}   double cdfErlang(long n, double b, double x)/* ============================================  * NOTE: use n >= 1, b > 0, and x > 0  * ============================================ */{   return (InGamma(n, x / b));}   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 (fabs(x - t) >= TINY);   return (x);}   static double pdfStandard(double x)/* ===================================  * NOTE: x can be any value  * =================================== */{   return (exp(- 0.5 * x * x) / SQRT2PI);}   static double cdfStandard(double x)/* ===================================  * NOTE: x can be any value  * =================================== */{