The R Project SVN R

Rev

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

Rev 59039 Rev 61535
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-2007 The R Core Team
 4 
 *  Copyright (C) 2000-2013 The R Core Team
5 
 *  Copyright (C) 2003	    The R Foundation
 5 
 *  Copyright (C) 2003-2013 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 44... Line 44...
44 
double qt(double p, double ndf, int lower_tail, int log_p)
 44 
double qt(double p, double ndf, int lower_tail, int log_p)
45 
{
 45 
{
46 
    const static double eps = 1.e-12;
 46 
    const static double eps = 1.e-12;
47 
 
 47 
 
48 
    double P, q;
 48 
    double P, q;
49 
    Rboolean neg;
 - 
 
50 
 
 49 
 
51 
#ifdef IEEE_754
 50 
#ifdef IEEE_754
52 
    if (ISNAN(p) || ISNAN(ndf))
 51 
    if (ISNAN(p) || ISNAN(ndf))
53 
	return p + ndf;
 52 
	return p + ndf;
54 
#endif
 53 
#endif
Line 60... Line 59...
60 
    if (ndf < 1) { /* based on qnt */
 59 
    if (ndf < 1) { /* based on qnt */
61 
	const static double accu = 1e-13;
 60 
	const static double accu = 1e-13;
62 
	const static double Eps = 1e-11; /* must be > accu */
 61 
	const static double Eps = 1e-11; /* must be > accu */
63 
 
 62 
 
64 
	double ux, lx, nx, pp;
 63 
	double ux, lx, nx, pp;
65
64 
 
66 
	int iter = 0;
 65 
	int iter = 0;
67 
 
 66 
 
68 
	p = R_DT_qIv(p);
 67 
	p = R_DT_qIv(p);
69 
 
 68 
 
70 
	/* Invert pt(.) :
 69 
	/* Invert pt(.) :
Line 100... Line 99...
100 
     */
 99 
     */
101 
    if (ndf > 1e20) return qnorm(p, 0., 1., lower_tail, log_p);
 100 
    if (ndf > 1e20) return qnorm(p, 0., 1., lower_tail, log_p);
102 
 
 101 
 
103 
    P = R_D_qIv(p); /* if exp(p) underflows, we fix below */
 102 
    P = R_D_qIv(p); /* if exp(p) underflows, we fix below */
104 
 
 103 
 
105 
    neg = (!lower_tail || P < 0.5) && (lower_tail || P > 0.5);
 104 
    Rboolean neg = (!lower_tail || P < 0.5) && (lower_tail || P > 0.5),
- 
 
 105 
	is_neg_lower = (lower_tail == neg); /* both TRUE or FALSE == !xor */
106 
    if(neg)
 106 
    if(neg)
107 
	P = 2 * (log_p ? (lower_tail ? P : -expm1(p)) : R_D_Lval(p));
 107 
	P = 2 * (log_p ? (lower_tail ? P : -expm1(p)) : R_D_Lval(p));
108 
    else
 108 
    else
109 
	P = 2 * (log_p ? (lower_tail ? -expm1(p) : P) : R_D_Cval(p));
 109 
	P = 2 * (log_p ? (lower_tail ? -expm1(p) : P) : R_D_Cval(p));
110 
    /* 0 <= P <= 1 ; P = 2*min(P', 1 - P')  in all cases */
 110 
    /* 0 <= P <= 1 ; P = 2*min(P', 1 - P')  in all cases */
111 
 
 111 
 
112 
/* Use this if(log_p) only : */
 - 
 
113 
#define P_is_exp_2p (lower_tail == neg) /* both TRUE or FALSE == !xor */
 - 
 
114 
 
 - 
 
115 
     if (fabs(ndf - 2) < eps) {	/* df ~= 2 */
 112 
     if (fabs(ndf - 2) < eps) {	/* df ~= 2 */
116 
	if(P > DBL_MIN) {
 113 
	if(P > DBL_MIN) {
117 
	    if(3* P < DBL_EPSILON) /* P ~= 0 */
 114 
	    if(3* P < DBL_EPSILON) /* P ~= 0 */
118 
		q = 1 / sqrt(P);
 115 
		q = 1 / sqrt(P);
119 
	    else if (P > 0.9)	   /* P ~= 1 */
 116 
	    else if (P > 0.9)	   /* P ~= 1 */
Line 121... Line 118...
121 
	    else /* eps/3 <= P <= 0.9 */
 118 
	    else /* eps/3 <= P <= 0.9 */
122 
		q = sqrt(2 / (P * (2 - P)) - 2);
 119 
		q = sqrt(2 / (P * (2 - P)) - 2);
123 
	}
 120 
	}
124 
	else { /* P << 1, q = 1/sqrt(P) = ... */
 121 
	else { /* P << 1, q = 1/sqrt(P) = ... */
125 
	    if(log_p)
 122 
	    if(log_p)
126 
		q = P_is_exp_2p ? exp(- p/2) / M_SQRT2 : 1/sqrt(-expm1(p));
 123 
		q = is_neg_lower ? exp(- p/2) / M_SQRT2 : 1/sqrt(-expm1(p));
127 
	    else
 124 
	    else
128 
		q = ML_POSINF;
 125 
		q = ML_POSINF;
129 
	}
 126 
	}
130 
    }
 127 
    }
131 
    else if (ndf < 1 + eps) { /* df ~= 1  (df < 1 excluded above): Cauchy */
 128 
    else if (ndf < 1 + eps) { /* df ~= 1  (df < 1 excluded above): Cauchy */
132 
	if(P > 0)
 129 
	if(P > 0)
133 
	    q = 1/tan(P * M_PI_2);/* == - tan((P+1) * M_PI_2) -- suffers for P ~= 0 */
 130 
	    q = 1/tan(P * M_PI_2);/* == - tan((P+1) * M_PI_2) -- suffers for P ~= 0 */
134 
 
 131 
 
135 
	else { /* P = 0, but maybe = 2*exp(p) ! */
 132 
	else { /* P = 0, but maybe = 2*exp(p) ! */
136 
	    if(log_p) /* 1/tan(e) ~ 1/e */
 133 
	    if(log_p) /* 1/tan(e) ~ 1/e */
137 
		q = P_is_exp_2p ? M_1_PI * exp(-p) : -1./(M_PI * expm1(p));
 134 
		q = is_neg_lower ? M_1_PI * exp(-p) : -1./(M_PI * expm1(p));
138 
	    else
 135 
	    else
139 
		q = ML_POSINF;
 136 
		q = ML_POSINF;
140 
	}
 137 
	}
141 
    }
 138 
    }
142 
    else {		/*-- usual case;  including, e.g.,  df = 1.1 */
 139 
    else {		/*-- usual case;  including, e.g.,  df = 1.1 */
Line 149... Line 146...
149 
	Rboolean P_ok1 = P > DBL_MIN || !log_p,  P_ok = P_ok1;
 146 
	Rboolean P_ok1 = P > DBL_MIN || !log_p,  P_ok = P_ok1;
150 
	if(P_ok1) {
 147 
	if(P_ok1) {
151 
	    y = pow(d * P, 2 / ndf);
 148 
	    y = pow(d * P, 2 / ndf);
152 
	    P_ok = (y >= DBL_EPSILON);
 149 
	    P_ok = (y >= DBL_EPSILON);
153 
	}
 150 
	}
154 
	if(!P_ok) { /* log_p && P very small */
 151 
	if(!P_ok) {// log.p && P very.small  ||  (d*P)^(2/df) =: y < eps_c
155 
	    log_P2 = P_is_exp_2p ? p : R_Log1_Exp(p); /* == log(P / 2) */
 152 
	    log_P2 = is_neg_lower ? R_D_log(p) : R_D_LExp(p); /* == log(P / 2) */
156 
	    x = (log(d) + M_LN2 + log_P2) / ndf;
 153 
	    x = (log(d) + M_LN2 + log_P2) / ndf;
157 
	    y = exp(2 * x);
 154 
	    y = exp(2 * x);
158 
	}
 155 
	}
159 
 
 156 
 
160 
	if ((ndf < 2.1 && P > 0.5) || y > 0.05 + a) { /* P > P0(df) */
 157 
	if ((ndf < 2.1 && P > 0.5) || y > 0.05 + a) { /* P > P0(df) */
Line 170... Line 167...
170 
	    c = (((0.05 * d * x - 5) * x - 7) * x - 2) * x + b + c;
 167 
	    c = (((0.05 * d * x - 5) * x - 7) * x - 2) * x + b + c;
171 
	    y = (((((0.4 * y + 6.3) * y + 36) * y + 94.5) / c
 168 
	    y = (((((0.4 * y + 6.3) * y + 36) * y + 94.5) / c
172 
		  - y - 3) / b + 1) * x;
 169 
		  - y - 3) / b + 1) * x;
173 
	    y = expm1(a * y * y);
 170 
	    y = expm1(a * y * y);
174 
	    q = sqrt(ndf * y);
 171 
	    q = sqrt(ndf * y);
175 
	} else { /* re-use 'y' from above */
 - 
 
176 
 
 - 
 
177 
	    if(!P_ok && x < - M_LN2 * DBL_MANT_DIG) {/* 0.5* log(DBL_EPSILON) */
 172 
	} else if(!P_ok && x < - M_LN2 * DBL_MANT_DIG) {/* 0.5* log(DBL_EPSILON) */
178 
		/* y above might have underflown */
 173 
	    /* y above might have underflown */
179 
		q = sqrt(ndf) * exp(-x);
 174 
	    q = sqrt(ndf) * exp(-x);
180 
	    }
 175 
	}
181 
	    else {
 176 
	else { /* re-use 'y' from above */
182 
		y = ((1 / (((ndf + 6) / (ndf * y) - 0.089 * d - 0.822)
 177 
	    y = ((1 / (((ndf + 6) / (ndf * y) - 0.089 * d - 0.822)
183 
			   * (ndf + 2) * 3) + 0.5 / (ndf + 4))
 178 
		       * (ndf + 2) * 3) + 0.5 / (ndf + 4))
184 
		     * y - 1) * (ndf + 1) / (ndf + 2) + 1 / y;
 179 
		 * y - 1) * (ndf + 1) / (ndf + 2) + 1 / y;
185 
		q = sqrt(ndf * y);
 180 
	    q = sqrt(ndf * y);
186 
	    }
 - 
 
187 
	}
 181 
	}
188 
 
 182 
 
189 
 
 183 
 
190 
	/* Now apply 2-term Taylor expansion improvement (1-term = Newton):
 184 
	/* Now apply 2-term Taylor expansion improvement (1-term = Newton):
191 
	 * as by Hill (1981) [ref.above] */
 185 
	 * as by Hill (1981) [ref.above] */