The R Project SVN R

Rev

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

Rev 68947 Rev 70234
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
3 
 *  Copyright (C) 1998 Ross Ihaka
 3 
 *  Copyright (C) 1998 Ross Ihaka
4 
 *  Copyright (C) 2000-2014 The R Core Team
 4 
 *  Copyright (C) 2000-2016 The R Core Team
5 
 *  Copyright (C) 2005 The R Foundation
 5 
 *  Copyright (C) 2005-2016 The R Foundation
6 
 *
 6 
 *
7 
 *  This program is free software; you can redistribute it and/or modify
 7 
 *  This program is free software; you can redistribute it and/or modify
8 
 *  it under the terms of the GNU General Public License as published by
 8 
 *  it under the terms of the GNU General Public License as published by
9 
 *  the Free Software Foundation; either version 2 of the License, or
 9 
 *  the Free Software Foundation; either version 2 of the License, or
10 
 *  (at your option) any later version.
 10 
 *  (at your option) any later version.
Line 45... Line 45...
45 
#include "dpq.h"
 45 
#include "dpq.h"
46 
 
 46 
 
47 
static double
 47 
static double
48 
do_search(double y, double *z, double p, double n, double pr, double incr)
 48 
do_search(double y, double *z, double p, double n, double pr, double incr)
49 
{
 49 
{
50 
    if(*z >= p) {
 - 
 
51 
			/* search to the left */
 50 
    if(*z >= p) {	/* search to the left */
52 
	for(;;) {
 51 
	for(;;) {
53 
	    if(y == 0 ||
 52 
	    if(y == 0 ||
54 
	       (*z = pnbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
 53 
	       (*z = pnbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
55 
		return y;
 54 
		return y;
56 
	    y = fmax2(0, y - incr);
 55 
	    y = fmax2(0, y - incr);
57 
	}
 56 
	}
58 
    }
 57 
    }
59 
    else {		/* search to the right */
 58 
    else {		/* search to the right */
60 
 
 - 
 
61 
	for(;;) {
 59 
	for(;;) {
62 
	    y = y + incr;
 60 
	    y = y + incr;
63 
	    if((*z = pnbinom(y, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) >= p)
 61 
	    if((*z = pnbinom(y, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) >= p)
64 
		return y;
 62 
		return y;
65 
	}
 63 
	}
Line 80... Line 78...
80 
       prob == size/(size+mu)
 78 
       prob == size/(size+mu)
81 
    */
 79 
    */
82 
    if (prob == 0 && size == 0) return 0;
 80 
    if (prob == 0 && size == 0) return 0;
83 
 
 81 
 
84 
    if (prob <= 0 || prob > 1 || size < 0) ML_ERR_return_NAN;
 82 
    if (prob <= 0 || prob > 1 || size < 0) ML_ERR_return_NAN;
85
83 
 
86 
    if (prob == 1 || size == 0) return 0;
 84 
    if (prob == 1 || size == 0) return 0;
87 
 
 85 
 
88 
    R_Q_P01_boundaries(p, 0, ML_POSINF);
 86 
    R_Q_P01_boundaries(p, 0, ML_POSINF);
89 
 
 87 
 
90 
    Q = 1.0 / prob;
 88 
    Q = 1.0 / prob;
Line 126... Line 124...
126 
    }
 124 
    }
127 
}
 125 
}
128 
 
 126 
 
129 
double qnbinom_mu(double p, double size, double mu, int lower_tail, int log_p)
 127 
double qnbinom_mu(double p, double size, double mu, int lower_tail, int log_p)
130 
{
 128 
{
- 
 
 129 
    if (size == ML_POSINF) // limit case: Poisson
- 
 
 130 
	return(qpois(p, mu, lower_tail, log_p));
131 
/* FIXME!  Implement properly!! (not losing accuracy for very large size (prob ~= 1)*/
 131 
/* FIXME!  Implement properly!! (not losing accuracy for very large size (prob ~= 1)*/
132 
    return qnbinom(p, size, /* prob = */ size/(size+mu), lower_tail, log_p);
 132 
    return qnbinom(p, size, /* prob = */ size/(size+mu), lower_tail, log_p);
133 
}
 133 
}