The R Project SVN R

Rev

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

Rev 18536 Rev 24748
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-2002 The R Development Core Team
 4 
 *  Copyright (C) 2000-2002 The R Development Core Team
- 
 
 5 
 *  Copyright (C) 2003	    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 22... Line 23...
22 
 *	The "Student" t distribution quantile function.
 23 
 *	The "Student" t distribution quantile function.
23 
 *
 24 
 *
24 
 *  NOTES
 25 
 *  NOTES
25 
 *
 26 
 *
26 
 *	This is a C translation of the Fortran routine given in:
 27 
 *	This is a C translation of the Fortran routine given in:
27 
 *	Algorithm 396: Student's t-quantiles by
 28 
 *	Hill, G.W (1970) "Algorithm 396: Student's t-quantiles"
28 
 *	G.W. Hill CACM 13(10), 619-620, October 1970
 29 
 *	CACM 13(10), 619-620.
- 
 
 30 
 *
- 
 
 31 
 *  ADDITIONS:
- 
 
 32 
 *	- lower_tail, log_p
- 
 
 33 
 *	- using	 expm1() : takes care of  Lozy (1979) "Remark on Algo.", TOMS
- 
 
 34 
 *	- Apply 2-term Taylor expansion as in
- 
 
 35 
 *	  Hill, G.W (1981) "Remark on Algo.396", ACM TOMS 7, 250-1
- 
 
 36 
 *	- Improve the formula decision for 1 < df < 2
29 
 */
 37 
 */
30 
 
 38 
 
31 
#include "nmath.h"
 39 
#include "nmath.h"
32 
#include "dpq.h"
 40 
#include "dpq.h"
33 
 
 41 
 
Line 59... Line 67...
59 
 
 67 
 
60 
    if((lower_tail && p_ > 0.5) || (!lower_tail && p_ < 0.5)) {
 68 
    if((lower_tail && p_ > 0.5) || (!lower_tail && p_ < 0.5)) {
61 
	neg = FALSE; P = 2 * R_D_Cval(p_);
 69 
	neg = FALSE; P = 2 * R_D_Cval(p_);
62 
    } else {
 70 
    } else {
63 
	neg = TRUE;  P = 2 * R_D_Lval(p_);
 71 
	neg = TRUE;  P = 2 * R_D_Lval(p_);
64 
    } /* 0 <= P <= 1  in all cases */
 72 
    } /* 0 <= P <= 1 ; P = 2*min(p_, 1 - p_)  in all cases */
65 
 
 73 
 
66 
    if (fabs(ndf - 2) < eps) {	/* df ~= 2 */
 74 
    if (fabs(ndf - 2) < eps) {	/* df ~= 2 */
67 
	if(P > 0)
 75 
	if(P > 0)
68 
	    q = sqrt(2 / (P * (2 - P)) - 2);
 76 
	    q = sqrt(2 / (P * (2 - P)) - 2);
69 
	else { /* P = 0, but maybe = exp(p) ! */
 77 
	else { /* P = 0, but maybe = exp(p) ! */
70 
	    if(log_p) q = M_SQRT2 * exp(- .5 * R_D_Lval(p));
 78 
	    if(log_p) q = M_SQRT2 * exp(- .5 * R_D_Lval(p));
71 
	    else q = ML_POSINF;
 79 
	    else q = ML_POSINF;
72 
	}
 80 
	}
73 
    }
 81 
    }
74 
    else if (ndf < 1 + eps) { /* df ~= 1  (df < 1 excluded above !) */
 82 
    else if (ndf < 1 + eps) { /* df ~= 1  (df < 1 excluded above): Cauchy */
75 
	if(P > 0)
 83 
	if(P > 0)
76 
	    q = - tan((P+1) * M_PI_2);
 84 
	    q = - tan((P+1) * M_PI_2);
77 
 
 85 
 
78 
	else { /* P = 0, but maybe p_ = exp(p) ! */
 86 
	else { /* P = 0, but maybe p_ = exp(p) ! */
79 
	    if(log_p) q = M_1_PI * exp(-R_D_Lval(p));/* cot(e) ~ 1/e */
 87 
	    if(log_p) q = M_1_PI * exp(-R_D_Lval(p));/* cot(e) ~ 1/e */
Line 85... Line 93...
85 
	b = 48 / (a * a);
 93 
	b = 48 / (a * a);
86 
	c = ((20700 * a / b - 98) * a - 16) * a + 96.36;
 94 
	c = ((20700 * a / b - 98) * a - 16) * a + 96.36;
87 
	d = ((94.5 / (b + c) - 3) / b + 1) * sqrt(a * M_PI_2) * ndf;
 95 
	d = ((94.5 / (b + c) - 3) / b + 1) * sqrt(a * M_PI_2) * ndf;
88 
	if(P > 0 || !log_p)
 96 
	if(P > 0 || !log_p)
89 
	    y = pow(d * P, 2 / ndf);
 97 
	    y = pow(d * P, 2 / ndf);
90 
	else /* P = 0 && log_p;  P = 2*exp(p*) */
 98 
	else /* P = 0 && log_p;	 P = 2*exp(p) */
91 
	    y = exp(2 / ndf * (log(d) + M_LN2 + R_D_Lval(p)));
 99 
	    y = exp(2 / ndf * (log(d) + M_LN2 + R_D_Lval(p)));
92 
 
 100 
 
93 
	if (y > 0.05 + a) {
 101 
	if ((ndf < 2.1 && P > 0.5) || y > 0.05 + a) { /* P > P0(df) */
94 
	    /* Asymptotic inverse expansion about normal */
 102 
	    /* Asymptotic inverse expansion about normal */
95 
	    if(P > 0 || !log_p)
 103 
	    if(P > 0 || !log_p)
96 
		x = qnorm(0.5 * P, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE);
 104 
		x = qnorm(0.5 * P, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE);
97 
	    else /* P = 0 && log_p;  P = 2*exp(p') */
 105 
	    else /* P = 0 && log_p;  P = 2*exp(p') */
98 
		x = qnorm( p,	   0., 1., lower_tail,	       /*log_p*/TRUE);
 106 
		x = qnorm( p,	   0., 1., lower_tail,	       /*log_p*/TRUE);
Line 108... Line 116...
108 
	    y = ((1 / (((ndf + 6) / (ndf * y) - 0.089 * d - 0.822)
 116 
	    y = ((1 / (((ndf + 6) / (ndf * y) - 0.089 * d - 0.822)
109 
		       * (ndf + 2) * 3) + 0.5 / (ndf + 4))
 117 
		       * (ndf + 2) * 3) + 0.5 / (ndf + 4))
110 
		 * y - 1) * (ndf + 1) / (ndf + 2) + 1 / y;
 118 
		 * y - 1) * (ndf + 1) / (ndf + 2) + 1 / y;
111 
	}
 119 
	}
112 
	q = sqrt(ndf * y);
 120 
	q = sqrt(ndf * y);
- 
 
 121 
 
- 
 
 122 
	/* Now apply 2-term Taylor expansion improvement (1-term = Newton):
- 
 
 123 
	 * as by Hill (1981) [ref.above] */
- 
 
 124 
 
- 
 
 125 
	/* FIXME: This is can be far from optimal when log_p = TRUE !
- 
 
 126 
	 *	  and probably also improvable when  lower_tail = FALSE */
- 
 
 127 
	x = (pt(q, ndf, /*lower_tail = */FALSE, /*log_p = */FALSE) - P/2) /
- 
 
 128 
	    dt(q, ndf, /* give_log = */FALSE);
- 
 
 129 
	/* Newton (=Taylor 1 term):
- 
 
 130 
	 *  q += x;
- 
 
 131 
	 * Taylor 2-term : */
- 
 
 132 
	q += x * (1. + x * q * (ndf + 1) / (2 * (q * q + ndf)));
113 
    }
 133 
    }
114 
    if(neg) q = -q;
 134 
    if(neg) q = -q;
115 
    return q;
 135 
    return q;
116 
}
 136 
}