The R Project SVN R

Rev

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

Rev 13688 Rev 27948
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  R : A Computer Language for Statistical Data Analysis
 2 
 *  R : A Computer Language for Statistical Data Analysis
3 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
 3 
 *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
4 
 *  Copyright (C) 2000, 2001 The R Development Core Team
 4 
 *  Copyright (C) 2000, 2001 The R Development Core Team
- 
 
 5 
 *  Copyright (C) 2004	     The R Foundation
5 
 *
 6 
 *
6 
 *  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
7 
 *  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
8 
 *  the Free Software Foundation; either version 2 of the License, or
 9 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 10 
 *  (at your option) any later version.
Line 11... Line 12...
11 
 *  This program is distributed in the hope that it will be useful,
 12 
 *  This program is distributed in the hope that it will be useful,
12 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 13 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
13 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 14 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 
 *  GNU General Public License for more details.
 15 
 *  GNU General Public License for more details.
15 
 *
 16 
 *
16 
 *  You should have received a copy of the GNU General Public License
 17 
 *  A copy of the GNU General Public License is available via WWW at
- 
 
 18 
 *  http://www.gnu.org/copyleft/gpl.html.  You can also obtain it by
17 
 *  along with this program; if not, write to the Free Software
 19 
 *  writing to the Free Software Foundation, Inc., 59 Temple Place,
18 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
 20 
 *  Suite 330, Boston, MA  02111-1307  USA.
19 
 */
 21 
 */
20 
 
 22 
 
21 
#include "nmath.h"
 23 
#include "nmath.h"
22 
#include "dpq.h"
 24 
#include "dpq.h"
23 
 
 25 
 
24 
double qnchisq(double p, double n, double lambda, int lower_tail, int log_p)
 26 
double qnchisq(double p, double n, double lambda, int lower_tail, int log_p)
25 
{
 27 
{
26 
    const double acu = 1e-12;
 28 
    const double acu = 1e-12;
27 
    const double Eps = 1e-6; /* must be > acu */
 29 
    const double Eps = 1e-6; /* must be > acu */
28 
 
 30 
 
29 
    double ux, lx, nx;
 31 
    double ux, lx, nx, pp;
30 
 
 32 
 
31 
#ifdef IEEE_754
 33 
#ifdef IEEE_754
32 
    if (ISNAN(p) || ISNAN(n) || ISNAN(lambda))
 34 
    if (ISNAN(p) || ISNAN(n) || ISNAN(lambda))
33 
	return p + n + lambda;
 35 
	return p + n + lambda;
34 
#endif
 36 
#endif
Line 37... Line 39...
37 
    n = floor(n + 0.5);
 39 
    n = floor(n + 0.5);
38 
    if (n < 1 || lambda < 0) ML_ERR_return_NAN;
 40 
    if (n < 1 || lambda < 0) ML_ERR_return_NAN;
39 
 
 41 
 
40 
    R_Q_P01_check(p);
 42 
    R_Q_P01_check(p);
41 
 
 43 
 
42 
    if (p == R_DT_0)
 44 
    if (p == R_DT_0) return 0;
43 
	return 0;
 - 
 
44 
 
 - 
 
45 
    /* Invert pnchisq(.) finding an upper and lower bound;
 - 
 
46 
       then interval halfing : */
 45 
    if (p == R_DT_1) return ML_POSINF;
47 
 
 46 
 
48 
    p = R_D_qIv(p);
 47 
    p = R_D_qIv(p);
- 
 
 48 
 
- 
 
 49 
    /* Invert pnchisq(.) :
- 
 
 50 
     * 1. finding an upper and lower bound */
49 
    if(lower_tail) {
 51 
    if(lower_tail) {
- 
 
 52 
	pp = p * (1 + Eps);/*not good when p ~= 1; caught via DBL_MAX */
- 
 
 53 
        for(ux = 1.;
50 
        for(ux = 1.; pnchisq_raw(ux, n, lambda, Eps, 128) < p * (1 + Eps);
54 
	    ux < DBL_MAX && pnchisq_raw(ux, n, lambda, Eps, 128) < pp;
51 
	    ux *= 2);
 55 
	    ux *= 2);
- 
 
 56 
	pp = p * (1 - Eps);
- 
 
 57 
        for(lx = fmin2(ux, DBL_MAX);
52 
        for(lx = ux; pnchisq_raw(lx, n, lambda, Eps, 128) > p * (1 - Eps);
58 
	    lx > DBL_MIN && pnchisq_raw(lx, n, lambda, Eps, 128) > pp;
53 
	    lx *= 0.5);
 59 
	    lx *= 0.5);
54 
    }
 60 
    }
55 
    else {
 61 
    else {
- 
 
 62 
	pp = (p > Eps) ? 1 + Eps : 1;
- 
 
 63 
        for(ux = 1.;
56 
        for(ux = 1.; pnchisq_raw(ux, n, lambda, Eps, 128) + p < 1 + Eps;
64 
	    ux < DBL_MAX && pnchisq_raw(ux, n, lambda, Eps, 128) + p < pp;
57 
	    ux *= 2);
 65 
	    ux *= 2);
- 
 
 66 
	pp = 1 - Eps;
- 
 
 67 
        for(lx = fmin2(ux, DBL_MAX);
58 
        for(lx = ux; pnchisq_raw(lx, n, lambda, Eps, 128) + p > 1 - Eps;
68 
	    lx > DBL_MIN && pnchisq_raw(lx, n, lambda, Eps, 128) + p > pp;
59 
	    lx *= 0.5);
 69 
	    lx *= 0.5);
60 
    }
 70 
    }
- 
 
 71 
 
61 
    p = R_D_Lval(p);
 72 
    p = R_D_Lval(p);
- 
 
 73 
 
- 
 
 74 
    /* 2. interval (lx,ux)  halving : */
62 
    do {
 75 
    do {
63 
	nx = 0.5 * (lx + ux);
 76 
	nx = 0.5 * (lx + ux);
64 
	if (pnchisq_raw(nx, n, lambda, acu, 1000) > p)
 77 
	if (pnchisq_raw(nx, n, lambda, acu, 1000) > p)
65 
	    ux = nx;
 78 
	    ux = nx;
66 
	else
 79 
	else