The R Project SVN R

Rev

Rev 46805 | Rev 64668 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
42303 ripley 1 
/*
2 
 *  AUTHOR
3 
 *    Catherine Loader, catherine@research.bell-labs.com.
4 
 *    October 23, 2000 and Feb, 2001.
5 
 *
46039 maechler 6 
 *    dnbinom_mu(): Martin Maechler, June 2008
7 
 *
42303 ripley 8 
 *  Merge in to R:
59037 ripley 9 
 *	Copyright (C) 2000--2008, The R Core Team
42303 ripley 10 
 *
11 
 *  This program is free software; you can redistribute it and/or modify
12 
 *  it under the terms of the GNU General Public License as published by
13 
 *  the Free Software Foundation; either version 2 of the License, or
14 
 *  (at your option) any later version.
15 
 *
16 
 *  This program is distributed in the hope that it will be useful,
17 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
18 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19 
 *  GNU General Public License for more details.
20 
 *
21 
 *  You should have received a copy of the GNU General Public License
22 
 *  along with this program; if not, a copy is available at
23 
 *  http://www.r-project.org/Licenses/
24 
 *
25 
 *
26 
 * DESCRIPTION
27 
 *
28 
 *   Computes the negative binomial distribution. For integer n,
29 
 *   this is probability of x failures before the nth success in a
30 
 *   sequence of Bernoulli trials. We do not enforce integer n, since
31 
 *   the distribution is well defined for non-integers,
32 
 *   and this can be useful for e.g. overdispersed discrete survival times.
33 
 */
34 
 
35 
#include "nmath.h"
36 
#include "dpq.h"
37 
 
44785 ripley 38 
double dnbinom(double x, double size, double prob, int give_log)
46039 maechler 39 
{
44785 ripley 40 
    double ans, p;
42303 ripley 41 
 
42 
#ifdef IEEE_754
44785 ripley 43 
    if (ISNAN(x) || ISNAN(size) || ISNAN(prob))
44 
        return x + size + prob;
42303 ripley 45 
#endif
46 
 
46039 maechler 47 
    if (prob <= 0 || prob > 1 || size < 0) ML_ERR_return_NAN;
42303 ripley 48 
    R_D_nonint_check(x);
49 
    if (x < 0 || !R_FINITE(x)) return R_D__0;
50 
    x = R_D_forceint(x);
51 
 
44785 ripley 52 
    ans = dbinom_raw(size, x+size, prob, 1-prob, give_log);
53 
    p = ((double)size)/(size+x);
54 
    return((give_log) ? log(p) + ans : p * ans);
42303 ripley 55 
}
46039 maechler 56 
 
57 
double dnbinom_mu(double x, double size, double mu, int give_log)
58 
{
46485 maechler 59 
    /* originally, just set  prob :=  size / (size + mu)  and called dbinom_raw(),
60 
     * but that suffers from cancellation when   mu << size  */
46039 maechler 61 
    double ans, p;
62 
 
63 
#ifdef IEEE_754
64 
    if (ISNAN(x) || ISNAN(size) || ISNAN(mu))
65 
        return x + size + mu;
66 
#endif
67 
 
68 
    if (mu < 0 || size < 0) ML_ERR_return_NAN;
69 
    R_D_nonint_check(x);
70 
    if (x < 0 || !R_FINITE(x)) return R_D__0;
71 
    x = R_D_forceint(x);
46485 maechler 72 
    if(x == 0)/* be accurate, both for n << mu, and n >> mu :*/
73 
	return R_D_exp(size * (size < mu ? log(size/(size+mu)) : log1p(- mu/(size+mu))));
46039 maechler 74 
    if(x < 1e-10 * size) { /* don't use dbinom_raw() but MM's formula: */
46805 maechler 75 
	/* FIXME --- 1e-8 shows problem; rather use algdiv() from ./toms708.c */
46039 maechler 76 
	return R_D_exp(x * log(size*mu / (size+mu)) - mu - lgamma(x+1) +
77 
		       log1p(x*(x-1)/(2*size)));
78 
    }
79 
    /* else: no unnecessary cancellation inside dbinom_raw, when
46485 maechler 80 
     * x_ = size and n_ = x+size are so close that n_ - x_ loses accuracy
46039 maechler 81 
     */
82 
    ans = dbinom_raw(size, x+size, size/(size+mu), mu/(size+mu), give_log);
83 
    p = ((double)size)/(size+x);
84 
    return((give_log) ? log(p) + ans : p * ans);
85 
}