lib/gsl-extras/binomial.c

   1 /* cdf/binomial.c
   2  *
   3  * Copyright (C) 2004 Free Software Foundation, Inc.
   4  *
   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 2 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  19  */
  20
  21 /*
  22  * Computes the cumulative distribution function for a binomial
  23  * random variable. For a binomial random variable X with n trials
  24  * and success probability p,
  25  *
  26  *          Pr( X <= k ) = Pr( Y >= p )
  27  *
  28  * where Y is a beta random variable with parameters k+1 and n-k.
  29  *
  30  * Reference:
  31  *
  32  * W. Feller, "An Introduction to Probability and Its
  33  * Applications," volume 1. Wiley, 1968. Exercise 45, page 173,
  34  * chapter 6.
  35  */
  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