lib/gsl-extras/geometric.c

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