confidence_intervals.c

来自「置信区间的相关程序」· C语言 代码 · 共 197 行

/* confidence_intervals.c
 *
 * Author:       Jon Hamkins
 * Date:         Tue Dec 11 09:45:02 PST 2001
 * Language:     ANSI C
 * Last Revised: Tue Dec 11 11:32:41 PST 2001
 *
 * Functions to compute confidence intervals.  Given x events occuring
 * in N tries, compute the 95% (or other) confidence interval, assuming
 * one of
 *
 * (1) a binomial distribution, useful when the events occur
 * independently at discrete times, or
 *
 * (2) a Poisson distribution, useful when the arrival times occur
 * independently (in continuous time), or as an approximation to a
 * binomial distribution when N->oo, x->constant.
 *
 * (3) a Gaussian distribution, useful as an approximation of binomial
 * distribution when N->oo and x/N->constant.  Accurate when the number
 * of trials is large.
 *
 * For example, the 95% confidence interval assuming a Gaussian
 * distribution for b errors in N simulated bits is the interval
 * (gaussian_low(b,N,0.025), gaussian_high(b,N,0.025)).
 *
 * Refs: 
 *
 * Stark & Woods "Probability, Random Processes, and Estimation
 * Theory for Engineers"
 *
 * http://members.aol.com/johnp71/confint.html
 *
 */

#include <math.h>

#define EPS (0.000001)
#define TOOBIG (1e30)

/* Compute probability of at least x1 successes, and no more than x2
 * successes, in N tries from a binomial distribution with population
 * proportion p. */
double
binomial_p(double N, double p, double x1, double x2)
{
   double q, k = 0., v = 1., s = 0., tot = 0.;

   q = p / (1. - p);
   while (k <= N) {
      tot += v;
      if ((k >= x1) && (k <= x2)) s += v;
      if (tot > TOOBIG) {
         s /= TOOBIG;
         tot /= TOOBIG;
         v /= TOOBIG;
      }
      k++;
      v *= q * (N + 1 - k) / k;
   }
   return(s/tot);
}

/* Perform binary search to find lower confidence interval, for binomial
 * distribution */
double 
binomial_low(double vx, double vN, double p)
{
   double vP, v, vsL = 0., vsH;

   vP = vx / vN;
   v = 0.5 * vP;
   vsH = vP;
   while ((vsH - vsL) > EPS) {
      if (binomial_p(vN, v, vx, vN) > p) {
         vsH = v;
         v = 0.5 * (vsL + v);
      }
      else {
         vsL = v;
         v = 0.5 * (v + vsH);
      }
   }
   return(v);
}

/* Perform binary search to find upper confidence interval, for binomial
 * distribution */
double
binomial_high(double vx, double vN, double p) 
{
   double vP, v, vsL, vsH = 1.;

   vP = vx / vN;
   v = 0.5 * (1. + vP);
   vsL = vP;
   while ((vsH - vsL) > EPS) {
      if (binomial_p(vN, v, 0., vx) < p) {
         vsH = v;
         v = 0.5 * (vsL + v);
      }
      else {
         vsL = v;
         v = 0.5 * (v + vsH);
      }
   }
   return(v);
}

/* Compute probability of at least x1 successes, and no more than x2
 * successes, from a Poisson distribution, given z successes observed. */
double
poisson_p(double z, double x1, double x2)
{
   double q=1, tot=0, s=0, k=0;

   while ((k <= z) || (q > tot * 0.0000000001)) {
      tot += q;
      if ((k >= x1) && (k <= x2)) s += q;
      if (tot > TOOBIG) {
         s /= TOOBIG;
         tot /= TOOBIG;
         q /= TOOBIG;
      }
      k++;
      q *= z / k;
   }
   return(s / tot);
}

/* Perform binary search to find lower confidence interval, for Poisson
 * distribution */
double
poisson_low(double vz, double p)
{
   double v = 0.5, dv = 0.5;

   while (dv > 0.0000001) {
      dv *= 0.5;
      if (poisson_p((1 + vz) * v / (1 - v), vz, 10000000000) > p)
         v -= dv;
      else v += dv;
   }
   return((1.+vz)*v/(1.-v));
}

/* Perform binary search to find higher confidence interval, for Poisson
 * distribution */
double
poisson_high(double vz, double p)
{
   double v=0.5, dv=0.5;

   while (dv > 0.0000001) {
    dv *= 0.5;
    if (poisson_p((1 + vz) * v / (1 - v), 0, vz) < p)
        v -= dv;
        else v += dv;
   }
   return((1 + vz) * v / (1 - v));
}

/* Find lower confidence interval, for Gaussian approximation to
 * binomial distribution.  (A binomial r.v. is a sum of Bournoulli
 * r.v.'s, and by central limit theorem, a sum with many terms
 * approaches a Gaussian r.v. with mean equal to the sum.)  */
double 
gaussian_low(double vx, double vN, double p)
{
   double low, med, high, mean, var;

   mean = vx;
   var = vx*(1-vx/vN);
   low = -1000.;
   med = 0.5 * vx;
   high = vx;

   while ((high - low) > EPS) {
      if (0.5+0.5*erf(M_SQRT1_2*(med-mean)/sqrt(var)) > p) {
         high = med;
         med = 0.5 * (low + med);
      }
      else {
         low = med;
         med = 0.5 * (med + high);
      }
   }
   return(med/vN);
}

/* Find upper confidence interval (by symmetry about mean) */
double
gaussian_high(double vx, double vN, double p) 
{
   return(2*vx/vN-gaussian_low(vx,vN,p));
}