lib/gsl-extras/binomial.c

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