htswanalysis/MACS/lib/gsl/gsl-1.11/cdf/binomial.c

   1 /* cdf/binomial.c
   2  *
   3  * Copyright (C) 2004 Jason H. Stover.
   4  * Copyright (C) 2004 Giulio Bottazzi
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  19  */
  20
  21 #include <config.h>
  22 #include <math.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_cdf.h>
  25 #include <gsl/gsl_sf_gamma.h>
  26
  27 #include "error.h"
  28
  29 /* Computes the cumulative distribution function for a binomial
  30    random variable. For a binomial random variable X with n trials
  31    and success probability p,
  32
  33            Pr( X <= k ) = Pr( Y >= p )
  34
  35    where Y is a beta random variable with parameters k+1 and n-k.
  36
  37    The binomial distribution has the form,
  38
  39    prob(k) =  n!/(k!(n-k)!) *  p^k (1-p)^(n-k) for k = 0, 1, ..., n
  40
  41    The cumulated distributions can be expressed in terms of normalized
  42    incomplete beta functions (see Abramowitz & Stegun eq. 26.5.26 and
  43    eq. 6.6.3).
  44
  45    Reference:
  46
  47    W. Feller, "An Introduction to Probability and Its
  48    Applications," volume 1. Wiley, 1968. Exercise 45, page 173,
  49    chapter 6.
  50  */
  51
  52 #include <config.h>
  53 #include <math.h>
  54 #include <gsl/gsl_math.h>
  55 #include <gsl/gsl_errno.h>
  56 #include <gsl/gsl_cdf.h>
  57
  58 double
  59 gsl_cdf_binomial_P (const unsigned int k, const double p, const unsigned int n)
  60 {
  61   double P;
  62   double a;
  63   double b;
  64
  65   if (p > 1.0 || p < 0.0)
  66     {
  67       CDF_ERROR ("p < 0 or p > 1", GSL_EDOM);
  68     }
  69
  70   if (k >= n)
  71     {
  72       P = 1.0;
  73     }
  74   else
  75     {
  76       a = (double) k + 1.0;
  77       b = (double) n - k;
  78       P = gsl_cdf_beta_Q (p, a, b);
  79     }
  80
  81   return P;
  82 }
  83
  84 double
  85 gsl_cdf_binomial_Q (const unsigned int k, const double p, const unsigned int n)
  86 {
  87   double Q;
  88   double a;
  89   double b;
  90
  91   if (p > 1.0 || p < 0.0)
  92     {
  93       CDF_ERROR ("p < 0 or p > 1", GSL_EDOM);
  94     }
  95
  96   if (k >= n)
  97     {
  98       Q = 0.0;
  99     }
 100   else
 101     {
 102       a = (double) k + 1.0;
 103       b = (double) n - k;
 104       Q = gsl_cdf_beta_P (p, a, b);
 105     }
 106
 107   return Q;
 108 }