| /* |
| * AUTHOR |
| * Catherine Loader, catherine@research.bell-labs.com. |
| * October 23, 2000 and Feb, 2001. |
| * |
| * dnbinom_mu(): Martin Maechler, June 2008 |
| * |
| * Merge in to R: |
| * Copyright (C) 2000--2016, The R Core Team |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| * |
| * |
| * DESCRIPTION |
| * |
| * Computes the negative binomial distribution. For integer n, |
| * this is probability of x failures before the nth success in a |
| * sequence of Bernoulli trials. We do not enforce integer n, since |
| * the distribution is well defined for non-integers, |
| * and this can be useful for e.g. overdispersed discrete survival times. |
| */ |
| |
| #include "nmath.h" |
| #include "dpq.h" |
| |
| double dnbinom(double x, double size, double prob, int give_log) |
| { |
| double ans, p; |
| |
| #ifdef IEEE_754 |
| if (ISNAN(x) || ISNAN(size) || ISNAN(prob)) |
| return x + size + prob; |
| #endif |
| |
| if (prob <= 0 || prob > 1 || size < 0) ML_ERR_return_NAN; |
| R_D_nonint_check(x); |
| if (x < 0 || !R_FINITE(x)) return R_D__0; |
| /* limiting case as size approaches zero is point mass at zero */ |
| if (x == 0 && size==0) return R_D__1; |
| x = R_forceint(x); |
| if(!R_FINITE(size)) size = DBL_MAX; |
| |
| ans = dbinom_raw(size, x+size, prob, 1-prob, give_log); |
| p = ((double)size)/(size+x); |
| return((give_log) ? log(p) + ans : p * ans); |
| } |
| |
| double dnbinom_mu(double x, double size, double mu, int give_log) |
| { |
| /* originally, just set prob := size / (size + mu) and called dbinom_raw(), |
| * but that suffers from cancellation when mu << size */ |
| |
| #ifdef IEEE_754 |
| if (ISNAN(x) || ISNAN(size) || ISNAN(mu)) |
| return x + size + mu; |
| #endif |
| |
| if (mu < 0 || size < 0) ML_ERR_return_NAN; |
| R_D_nonint_check(x); |
| if (x < 0 || !R_FINITE(x)) return R_D__0; |
| |
| /* limiting case as size approaches zero is point mass at zero, |
| * even if mu is kept constant. limit distribution does not |
| * have mean mu, though. |
| */ |
| if (x == 0 && size == 0) return R_D__1; |
| x = R_forceint(x); |
| if(!R_FINITE(size)) // limit case: Poisson |
| return(dpois_raw(x, mu, give_log)); |
| |
| if(x == 0)/* be accurate, both for n << mu, and n >> mu :*/ |
| return R_D_exp(size * (size < mu ? log(size/(size+mu)) : log1p(- mu/(size+mu)))); |
| if(x < 1e-10 * size) { /* don't use dbinom_raw() but MM's formula: */ |
| /* FIXME --- 1e-8 shows problem; rather use algdiv() from ./toms708.c */ |
| double p = (size < mu ? log(size/(1 + size/mu)) : log(mu / (1 + mu/size))); |
| return R_D_exp(x * p - mu - lgamma(x+1) + |
| log1p(x*(x-1)/(2*size))); |
| } else { |
| /* no unnecessary cancellation inside dbinom_raw, when |
| * x_ = size and n_ = x+size are so close that n_ - x_ loses accuracy */ |
| double p = ((double)size)/(size+x), |
| ans = dbinom_raw(size, x+size, size/(size+mu), mu/(size+mu), give_log); |
| return((give_log) ? log(p) + ans : p * ans); |
| } |
| } |