The R Project SVN R

Rev

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

Rev Author Line No. Line
574 ihaka 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
3 
 *  Copyright (C) 1998 Ross Ihaka
59039 ripley 4 
 *  Copyright (C) 2000-8 The R Core Team
574 ihaka 5 
 *
6 
 *  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
8 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
10 
 *
11 
 *  This program is distributed in the hope that it will be useful,
12 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
13 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 
 *  GNU General Public License for more details.
15 
 *
16 
 *  You should have received a copy of the GNU General Public License
42302 ripley 17 
 *  along with this program; if not, a copy is available at
18 
 *  http://www.r-project.org/Licenses/
574 ihaka 19 
 *
20 
 *  DESCRIPTION
21 
 *
7671 maechler 22 
 *	The distribution function of the negative binomial distribution.
574 ihaka 23 
 *
24 
 *  NOTES
25 
 *
7671 maechler 26 
 *	x = the number of failures before the n-th success
574 ihaka 27 
 */
28 
 
8431 ripley 29 
#include "nmath.h"
8074 maechler 30 
#include "dpq.h"
574 ihaka 31 
 
44785 ripley 32 
double pnbinom(double x, double size, double prob, int lower_tail, int log_p)
574 ihaka 33 
{
580 ihaka 34 
#ifdef IEEE_754
44785 ripley 35 
    if (ISNAN(x) || ISNAN(size) || ISNAN(prob))
36 
	return x + size + prob;
37 
    if(!R_FINITE(size) || !R_FINITE(prob))	ML_ERR_return_NAN;
580 ihaka 38 
#endif
44785 ripley 39 
    if (size <= 0 || prob <= 0 || prob > 1)	ML_ERR_return_NAN;
7671 maechler 40 
 
41 
    if (x < 0) return R_DT_0;
7867 ripley 42 
    if (!R_FINITE(x)) return R_DT_1;
37563 maechler 43 
    x = floor(x + 1e-7);
44785 ripley 44 
    return pbeta(prob, size, x + 1, lower_tail, log_p);
574 ihaka 45 
}
46039 maechler 46 
 
47 
double pnbinom_mu(double x, double size, double mu, int lower_tail, int log_p)
48 
{
49 
#ifdef IEEE_754
50 
    if (ISNAN(x) || ISNAN(size) || ISNAN(mu))
51 
	return x + size + mu;
52 
    if(!R_FINITE(size) || !R_FINITE(mu))	ML_ERR_return_NAN;
53 
#endif
54 
    if (size <= 0 || mu < 0)	ML_ERR_return_NAN;
55 
 
56 
    if (x < 0) return R_DT_0;
57 
    if (!R_FINITE(x)) return R_DT_1;
58 
    x = floor(x + 1e-7);
59 
    /* return
60 
     * pbeta(pr, size, x + 1, lower_tail, log_p);  pr = size/(size + mu), 1-pr = mu/(size+mu)
61 
     *
62 
     *= pbeta_raw(pr, size, x + 1, lower_tail, log_p)
63 
     *            x.  pin   qin
64 
     *=  bratio (pin,  qin, x., 1-x., &w, &wc, &ierr, log_p),  and return w or wc ..
65 
     *=  bratio (size, x+1, pr, 1-pr, &w, &wc, &ierr, log_p) */
66 
    {
67 
	int ierr;
68 
	double w, wc;
69 
	bratio(size, x+1, size/(size+mu), mu/(size+mu), &w, &wc, &ierr, log_p);
70 
	if(ierr)
71 
	    MATHLIB_WARNING(_("pnbinom_mu() -> bratio() gave error code %d"), ierr);
72 
	return lower_tail ? w : wc;
73 
    }
74 
}