The R Project SVN R

Rev

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

Rev 46039 Rev 46485
Line 54... Line 54...
54 
    return((give_log) ? log(p) + ans : p * ans);
 54 
    return((give_log) ? log(p) + ans : p * ans);
55 
}
 55 
}
56 
 
 56 
 
57 
double dnbinom_mu(double x, double size, double mu, int give_log)
 57 
double dnbinom_mu(double x, double size, double mu, int give_log)
58 
{
 58 
{
- 
 
 59 
    /* originally, just set  prob :=  size / (size + mu)  and called dbinom_raw(),
- 
 
 60 
     * but that suffers from cancellation when   mu << size  */
59 
    double ans, p;
 61 
    double ans, p;
60 
 
 62 
 
61 
#ifdef IEEE_754
 63 
#ifdef IEEE_754
62 
    if (ISNAN(x) || ISNAN(size) || ISNAN(mu))
 64 
    if (ISNAN(x) || ISNAN(size) || ISNAN(mu))
63 
        return x + size + mu;
 65 
        return x + size + mu;
Line 65... Line 67...
65 
 
 67 
 
66 
    if (mu < 0 || size < 0) ML_ERR_return_NAN;
 68 
    if (mu < 0 || size < 0) ML_ERR_return_NAN;
67 
    R_D_nonint_check(x);
 69 
    R_D_nonint_check(x);
68 
    if (x < 0 || !R_FINITE(x)) return R_D__0;
 70 
    if (x < 0 || !R_FINITE(x)) return R_D__0;
69 
    x = R_D_forceint(x);
 71 
    x = R_D_forceint(x);
70 
    if(x == 0)
 72 
    if(x == 0)/* be accurate, both for n << mu, and n >> mu :*/
71 
	return R_D_exp(size * log1p(- mu/(size+mu)));
 73 
	return R_D_exp(size * (size < mu ? log(size/(size+mu)) : log1p(- mu/(size+mu))));
72 
    if(x < 1e-10 * size) { /* don't use dbinom_raw() but MM's formula: */
 74 
    if(x < 1e-10 * size) { /* don't use dbinom_raw() but MM's formula: */
73 
	return R_D_exp(x * log(size*mu / (size+mu)) - mu - lgamma(x+1) +
 75 
	return R_D_exp(x * log(size*mu / (size+mu)) - mu - lgamma(x+1) +
74 
		       log1p(x*(x-1)/(2*size)));
 76 
		       log1p(x*(x-1)/(2*size)));
75 
    }
 77 
    }
76 
    /* else: no unnecessary cancellation inside dbinom_raw, when
 78 
    /* else: no unnecessary cancellation inside dbinom_raw, when
77 
     * x_ = size and n_ = x+size are so cloase that n_ - x_ loses accuracy
 79 
     * x_ = size and n_ = x+size are so close that n_ - x_ loses accuracy
78 
     */
 80 
     */
79 
    ans = dbinom_raw(size, x+size, size/(size+mu), mu/(size+mu), give_log);
 81 
    ans = dbinom_raw(size, x+size, size/(size+mu), mu/(size+mu), give_log);
80 
    p = ((double)size)/(size+x);
 82 
    p = ((double)size)/(size+x);
81 
    return((give_log) ? log(p) + ans : p * ans);
 83 
    return((give_log) ? log(p) + ans : p * ans);
82 
}
 84 
}