The R Project SVN R

Rev

Rev 11140 | Rev 36820 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 11140 Rev 11723
Line 1... Line 1...
1 
/*
 1 
/*
- 
 
 2 
 * AUTHOR
2 
 *  Mathlib : A C Library of Special Functions
 3 
 *   Catherine Loader, catherine@research.bell-labs.com.
- 
 
 4 
 *   October 23, 2000.
- 
 
 5 
 *
3 
 *  Copyright (C) 1998 Ross Ihaka
 6 
 *  Merge in to R:
4 
 *  Copyright (C) 2000 The R Development Core Team
 7 
 *	Copyright (C) 2000, The R Core Development Team
5 
 *
 8 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 9 
 *  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
 10 
 *  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
 11 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 12 
 *  (at your option) any later version.
Line 13... Line 16...
13 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 16 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 
 *  GNU General Public License for more details.
 17 
 *  GNU General Public License for more details.
15 
 *
 18 
 *
16 
 *  You should have received a copy of the GNU General Public License
 19 
 *  You should have received a copy of the GNU General Public License
17 
 *  along with this program; if not, write to the Free Software
 20 
 *  along with this program; if not, write to the Free Software
18 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 21 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
19 
 *
 - 
 
20 
 *  DESCRIPTION
 - 
 
21 
 *
 22 
 *
22 
 *    The density of the binomial distribution.
 - 
 
23 
 *
 23 
 *
24 
 * Using the new algorithm of C.Loader(1999) :
 24 
 * DESCRIPTION
25 
 *
 25 
 *
26 
 * This code provides functions for computing binomial probabilities, and
 26 
 *   To compute the binomial probability, call dbinom(x,n,p).
27 
 * attempts to be accurate for a full range of parameter values
 - 
 
28 
 * (standard algorithms are often inaccurate with large parameters).
 27 
 *   This checks for argument validity, and calls dbinom_raw().
29 
 *
 28 
 *
- 
 
 29 
 *   dbinom_raw() does the actual computation; note this is called by
30 
 * Merge in to R:
 30 
 *   other functions in addition to dbinom()).
- 
 
 31 
 *     (1) dbinom_raw() has both p and q arguments, when one may be represented
- 
 
 32 
 *         more accurately than the other (in particular, in df()).
- 
 
 33 
 *     (2) dbinom_raw() does NOT check that inputs x and n are integers. This
31 
 * Copyright (C) 2000, The R Core Development Team
 34 
 *         should be done in the calling function, where necessary.
32 
 *
 - 
 
33 
 * NOTE: Loader's original code is now split (and merged into R's extras)
 35 
 *     (3) Also does not check for 0 <= p <= 1 and 0 <= q <= 1 or NaN's.
34 
 *       into  ./dbinom.c, ./dpois.c and ./stirlerr.c
 36 
 *         Do this in the calling function.
35 
 */
 37 
 */
36 
 
 38 
 
37 
#include "nmath.h"
 39 
#include "nmath.h"
38 
#include "dpq.h"
 40 
#include "dpq.h"
39 
 
 41 
 
- 
 
 42 
double dbinom_raw(double x, double n, double p, double q, int give_log)
- 
 
 43 
{
- 
 
 44 
    double f, lc;
- 
 
 45 
 
- 
 
 46 
    if (p == 0) return((x == 0) ? R_D__1 : R_D__0);
- 
 
 47 
    if (q == 0) return((x == n) ? R_D__1 : R_D__0);
- 
 
 48 
 
- 
 
 49 
    if (x == 0) {
- 
 
 50 
	if(n == 0) return R_D__1;
- 
 
 51 
	lc = (p < 0.1) ? -bd0(n,n*q) - n*p : n*log(q);
- 
 
 52 
	return( R_D_exp(lc) );
- 
 
 53 
    }
- 
 
 54 
    if (x == n) {
- 
 
 55 
	lc = (q < 0.1) ? -bd0(n,n*p) - n*q : n*log(p);
- 
 
 56 
	return( R_D_exp(lc) );
- 
 
 57 
    }
- 
 
 58 
    if (x < 0 || x > n) return( R_D__0 );
- 
 
 59 
 
- 
 
 60 
    lc = stirlerr(n) - stirlerr(x) - stirlerr(n-x) - bd0(x,n*p) - bd0(n-x,n*q);
- 
 
 61 
    f = (M_2PI*x*(n-x))/n;
- 
 
 62 
 
- 
 
 63 
    return R_D_fexp(f,lc);
- 
 
 64 
}
- 
 
 65 
 
40 
double dbinom(double x, double n, double p, int give_log)
 66 
double dbinom(double x, double n, double p, int give_log)
41 
{
 67 
{
42 
    double lc;
 - 
 
43 
#ifdef IEEE_754
 68 
#ifdef IEEE_754
44 
    /* NaNs propagated correctly */
 69 
    /* NaNs propagated correctly */
45 
    if (ISNAN(x) || ISNAN(n) || ISNAN(p)) return x + n + p;
 70 
    if (ISNAN(x) || ISNAN(n) || ISNAN(p)) return x + n + p;
46 
#endif
 71 
#endif
47 
    n = floor(n + 0.5);
 - 
 
48 
    if(n < 0 || p < 0 || p > 1) ML_ERR_return_NAN;
 - 
 
49 
 
 72 
 
50 
    if(fabs(x - floor(x + 0.5)) > 1e-7) {
 73 
  if (p < 0 || p > 1 || R_D_notnnegint(n)) ML_ERR_return_NAN;
51 
	MATHLIB_WARNING("non-integer x = %f", x);
 - 
 
52 
	return R_D__0;
 74 
  R_D_nonint_check(x);
53 
    }
 75
54 
    if (x < 0 || x > n)
 - 
 
55 
	return R_D__0;
 - 
 
56 
    if (x == 0)
 - 
 
57 
	return (n == 0 || p == 0) ? R_D__1
 - 
 
58 
	    : (p == 1) ? R_D__0 : R_D_exp(n*log1p(-p));
 - 
 
59 
    /* x > 0 : */
 - 
 
60 
    if (p == 0 || p == 1)
 - 
 
61 
	return (x == n && p == 1) ? R_D__1 : R_D__0;
 - 
 
62 
    /* 0 < p < 1 : */
 - 
 
63 
    if (x == n)
 76 
  n = R_D_forceint(n);
64 
	return give_log ? n*log(p) : pow(p,n);/* or R_pow_di() {w/o checks}*/
 - 
 
65 
    /* else (0 < x < n) : */
 - 
 
66 
#ifndef OLD_dbinom
 - 
 
67 
    lc = stirlerr(n) - stirlerr(x) - stirlerr(n-x)
 - 
 
68 
	- bd0(x, n*p)
 - 
 
69 
	- bd0(n-x, n*(1.-p));
 - 
 
70 
    if (give_log)
 77 
  x = R_D_forceint(x);
71 
	return lc - M_LN_SQRT_2PI + .5*log(n/(x*(n-x)));
 - 
 
72 
    else
 - 
 
73 
	return exp(lc) * sqrt(n/(2*M_PI*x*(n-x)));
 - 
 
74 
#else
 78 
 
75 
    return R_D_exp(lfastchoose(n, x) + log(p) * x + (n - x) * log1p(-p));
 79 
  return dbinom_raw(x,n,p,1-p,give_log);
76 
#endif
 - 
 
77 
}
 80 
}