lib/gsl-extras/binomial.c

   1 /* cdf/binomial.c
   2  *
   3  * Copyright (C) 2004 Free Software Foundation, Inc.
   4  * Written by Jason H. Stover.
   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 <config.h>
  37 #include <math.h>
  38 #include <gsl/gsl_math.h>
  39 #include <gsl/gsl_errno.h>
  40 #include <gsl/gsl_cdf.h>
  41 #include "gsl-extras.h"
  42
  43 double
  44 gslextras_cdf_binomial_P(const long k, const long n, const double p)
  45 {
  46   double P;
  47   double a;
  48   double b;
  49
  50   if(p > 1.0 || p < 0.0)
  51     {
  52       GSLEXTRAS_CDF_ERROR("p < 0 or p > 1",GSL_EDOM);
  53     }
  54   if ( k >= n )
  55     {
  56       P = 1.0;
  57     }
  58   else if (k < 0)
  59     {
  60       P = 0.0;
  61     }
  62   else
  63     {
  64       a = (double) k+1;
  65       b = (double) n - k;
  66       P = gsl_cdf_beta_Q( p, a, b);
  67     }
  68
  69   return P;
  70 }
  71 double
  72 gslextras_cdf_binomial_Q(const long k, const long n, const double q)
  73 {
  74   double P;
  75   double a;
  76   double b;
  77
  78   if(q > 1.0 || q < 0.0)
  79     {
  80       GSLEXTRAS_CDF_ERROR("p < 0 or p > 1",GSL_EDOM);
  81     }
  82   if( k >= n )
  83     {
  84       P = 0.0;
  85     }
  86   else if ( k < 0 )
  87     {
  88       P = 1.0;
  89     }
  90   else
  91     {
  92       a = (double) k+1;
  93       b = (double) n - k;
  94       P = gsl_cdf_beta_P(q, a, b);
  95     }
  96
  97   return P;
  98 }
  99