The R Project SVN R

Rev

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

Rev 10745 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 15... Line 18...
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 
 *
 22 
 *
20 
 *  DESCRIPTION
 - 
 
21 
 *
 - 
 
22 
 *    The density function of the Poisson distribution.
 - 
 
23 
 *
 23 
 *
24 
 * Using the new algorithm of Clive Loader(1999) :
 24 
 * DESCRIPTION
25 
 *
 25 
 *
26 
 * The author of this software is Clive Loader, clive@bell-labs.com.
 - 
 
27 
 * Copyright (c) 1999-2000 Lucent Technologies, Bell Laboratories.
 - 
 
28 
 * Permission to use, copy, modify, and distribute this software for any
 - 
 
29 
 * purpose without fee is hereby granted, with the exceptions noted below,
 - 
 
30 
 * and provided that this entire notice is included in all copies of any
 26 
 *    dpois() checks argument validity and calls dpois_raw().
31 
 * software which is or includes a copy or modification of this software
 - 
 
32 
 * and in all copies of the supporting documentation for such software.
 - 
 
33 
 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
 - 
 
34 
 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT TECHNOLOGIES
 - 
 
35 
 * MAKE ANY REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE
 - 
 
36 
 * MERCHANTABILITY OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
 - 
 
37 
 *
 - 
 
38 
 * This code provides a function for computing Poisson probabilities, and
 - 
 
39 
 * attempts to be accurate for a full range of parameter values
 - 
 
40 
 * (standard algorithms are often inaccurate with large parameters).
 - 
 
41 
 *
 27 
 *
- 
 
 28 
 *    dpois_raw() computes the Poisson probability  lb^x exp(-lb) / x!.
42 
 * Merge in to R:
 29 
 *      This does not check that x is an integer, since dgamma() may
- 
 
 30 
 *      call this with a fractional x argument. Any necessary argument
43 
 * Copyright (C) 2000, The R Core Development Team
 31 
 *      checks should be done in the calling function.
44 
 *
 32 
 *
45 
 * NOTE: Loader's original code is now split (and merged into R's extras)
 - 
 
46 
 *       into  ./dbinom.c, ./dpois.c and ./stirlerr.c
 - 
 
47 
 */
 33 
 */
48 
 
 34 
 
49 
#include "nmath.h"
 35 
#include "nmath.h"
50 
#include "dpq.h"
 36 
#include "dpq.h"
51 
 
 37 
 
52 
double dpois(double x, double lambda, int give_log)
 38 
double dpois_raw(double x, double lambda, int give_log)
53 
{
 39 
{
- 
 
 40 
    if (lambda == 0) return( (x == 0) ? R_D__1 : R_D__0 );
- 
 
 41 
    if (x == 0) return( R_D_exp(-lambda) );
- 
 
 42 
    if (x < 0)  return( R_D__0 );
54 
 
 43 
 
- 
 
 44 
    return(R_D_fexp( M_2PI*x, -stirlerr(x)-bd0(x,lambda) ));
- 
 
 45 
}
- 
 
 46 
 
- 
 
 47 
double dpois(double x, double lambda, int give_log)
- 
 
 48 
{
55 
#ifdef IEEE_754
 49 
#ifdef IEEE_754
56 
    if(ISNAN(x) || ISNAN(lambda))
 50 
    if(ISNAN(x) || ISNAN(lambda))
57 
	return x + lambda;
 51 
        return x + lambda;
58 
#endif
 - 
 
59 
    if(fabs(x - floor(x + 0.5)) > 1e-7) {
 - 
 
60 
	MATHLIB_WARNING("non-integer x = %f", x);
 - 
 
61 
	return R_D__0;
 - 
 
62 
    }
 - 
 
63 
    if(lambda < 0.0) ML_ERR_return_NAN;
 - 
 
64 
 
 - 
 
65 
    if (x < 0 || !R_FINITE(x))
 - 
 
66 
	return R_D__0;
 - 
 
67 
    if (lambda == 0)
 - 
 
68 
	return (x > 0) ? R_D__0 : R_D__1 ;
 - 
 
69 
    if (x == 0)
 - 
 
70 
	return R_D_exp(-lambda);
 - 
 
71 
    /* else, normal case : */
 - 
 
72 
#ifndef OLD_dpois
 - 
 
73 
    if (give_log)
 - 
 
74 
	return - (stirlerr(x) + bd0(x,lambda) + .5 * log(x) + M_LN_SQRT_2PI);
 - 
 
75 
    else
 - 
 
76 
	return exp(-stirlerr(x) - bd0(x,lambda)) / sqrt(2*M_PI*x);
 - 
 
77 
#else
 - 
 
78 
    return R_D_exp(x * log(lambda) - lambda - lgammafn(x + 1));
 - 
 
79 
#endif
 52 
#endif
80 
}
 - 
 
81 
 
 53 
 
- 
 
 54 
    if (lambda < 0) ML_ERR_return_NAN;
- 
 
 55 
    R_D_nonint_check(x);
- 
 
 56 
    if (x < 0 || !R_FINITE(x)) return R_D__0;
- 
 
 57 
    x = R_D_forceint(x);
- 
 
 58 
 
- 
 
 59 
    return( dpois_raw(x,lambda,give_log) );
- 
 
 60 
}