lib/gsl-extras/geometric.c

   1 /* cdf/geometric.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  * Pr(X <= n) for a negative binomial random variable X, i.e.,
  23  * the probability of n or fewer failuers before success k.
  24  */
  25
  26 #include <math.h>
  27 #include <gsl/gsl_math.h>
  28 #include <gsl/gsl_errno.h>
  29 #include <gsl/gsl_sf.h>
  30 #include <gsl/gsl_cdf.h>
  31 #include "gsl-extras.h"
  32
  33 /*
  34  * Pr (X <= n), i.e., the probability of n or fewer
  35  * failures until the first success.
  36  */
  37 double
  38 gslextras_cdf_geometric_P (const long n, const double p)
  39 {
  40   double P;
  41   double a;
  42   int i;
  43   int m;
  44   double sign = 1.0;
  45   double term;
  46   double q;
  47
  48   if(p > 1.0 || p < 0.0)
  49     {
  50       GSLEXTRAS_CDF_ERROR("p < 0 or p > 1",GSL_EDOM);
  51     }
  52   if ( n < 0 )
  53     {
  54       return 0.0;
  55     }
  56   q = 1.0 - p;
  57   a = (double) n+1;
  58   if( p < GSL_DBL_EPSILON )
  59     {
  60       /*
  61        * 1.0 - pow(q,a) will overflow, so use
  62        * a Taylor series.
  63        */
  64       i = 2;
  65       m = n+1;
  66       term = exp(log(a) + log(p));
  67       P = term;
  68       while ( term > GSL_DBL_MIN && i < m)
  69         {
  70           term = exp (sign * gsl_sf_lnchoose(m,i) + i * log(p));
  71           P += term;
  72           i++;
  73           sign = -sign;
  74         }
  75     }
  76   else
  77     {
  78       P = 1.0 - pow ( q, a);
  79     }
  80   return P;
  81 }
  82 double
  83 gslextras_cdf_geometric_Q ( const long n, const double p)
  84 {
  85   double P;
  86   double q;
  87   double a;
  88
  89   if(p > 1.0 || p < 0.0)
  90     {
  91       GSLEXTRAS_CDF_ERROR("p < 0 or p > 1",GSL_EDOM);
  92     }
  93   if ( n < 0 )
  94     {
  95       P = 1.0;
  96     }
  97   else
  98     {
  99       a = (double) n+1;
 100       q = 1.0 - p;
 101       P = pow(q, a);
 102     }
 103
 104   return P;
 105 }