The R Project SVN R

Rev

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

Rev 11248 Rev 24735
Line 71... Line 71...
71 
    REprintf("qnorm(p=%10.7g, m=%g, s=%g, l.t.= %d, log= %d): q = %g\n",
 71 
    REprintf("qnorm(p=%10.7g, m=%g, s=%g, l.t.= %d, log= %d): q = %g\n",
72 
	     p,mu,sigma, lower_tail, log_p, q);
 72 
	     p,mu,sigma, lower_tail, log_p, q);
73 
#endif
 73 
#endif
74 
 
 74 
 
75 
 
 75 
 
76 
#ifdef OLD_qnorm
 - 
 
77 
    /* --- use  AS 111 --- */
 - 
 
78 
    if (fabs(q) <= 0.42) {
 - 
 
79 
 
 - 
 
80 
	/* 0.08 <= p <= 0.92 */
 - 
 
81 
 
 - 
 
82 
	r = q * q;
 - 
 
83 
	val = q * (((-25.44106049637 * r + 41.39119773534) * r
 - 
 
84 
		    - 18.61500062529) * r + 2.50662823884)
 - 
 
85 
	    / ((((3.13082909833 * r - 21.06224101826) * r
 - 
 
86 
		 + 23.08336743743) * r + -8.47351093090) * r + 1.0);
 - 
 
87 
    }
 - 
 
88 
    else {
 - 
 
89 
 
 - 
 
90 
	/* p < 0.08 or p > 0.92, set r = min(p, 1 - p) */
 - 
 
91 
 
 - 
 
92 
	if (q > 0)
 - 
 
93 
	    r = R_DT_CIv(p);/* 1-p */
 - 
 
94 
	else
 - 
 
95 
	    r = p_;/* = R_DT_Iv(p) ^=  p */
 - 
 
96 
#ifdef DEBUG_qnorm
 - 
 
97 
	REprintf("\t 'middle p': r = %7g\n", r);
 - 
 
98 
#endif
 - 
 
99 
 
 - 
 
100 
	if(r > DBL_EPSILON) {
 - 
 
101 
	    r = sqrt(- ((log_p &&
 - 
 
102 
			 ((lower_tail && q <= 0) || (!lower_tail && q > 0))) ?
 - 
 
103 
			p : /* else */ log(r)));
 - 
 
104 
#ifdef DEBUG_qnorm
 - 
 
105 
	    REprintf("\t new r = %7g ( =? sqrt(- log(r)) )\n", r);
 - 
 
106 
#endif
 - 
 
107 
	    val = (((2.32121276858 * r + 4.85014127135) * r
 - 
 
108 
		    - 2.29796479134) * r - 2.78718931138)
 - 
 
109 
		/ ((1.63706781897 * r + 3.54388924762) * r + 1.0);
 - 
 
110 
	    if (q < 0)
 - 
 
111 
		val = -val;
 - 
 
112 
	}
 - 
 
113 
	else if(r >= DBL_MIN) { /* r = p <= eps : Use Wichura */
 - 
 
114 
	    val = -2 * (log_p ? R_D_Lval(p) : log(R_D_Lval(p)));
 - 
 
115 
	    r = log(2 * M_PI * val);
 - 
 
116 
#ifdef DEBUG_qnorm
 - 
 
117 
	    REprintf("\t DBL_MIN <= r <= DBL_EPS: val = %g, new r = %g\n",
 - 
 
118 
		     val, r);
 - 
 
119 
#endif
 - 
 
120 
	    p = val * val;
 - 
 
121 
	    r = r/val + (2 - r)/p + (-14 + 6 * r - r * r)/(2 * p * val);
 - 
 
122 
	    val = sqrt(val * (1 - r));
 - 
 
123 
	    if(q < 0.0)
 - 
 
124 
		val = -val;
 - 
 
125 
	    return mu + sigma * val;
 - 
 
126 
	}
 - 
 
127 
	else {
 - 
 
128 
#ifdef DEBUG_qnorm
 - 
 
129 
	    REprintf("\t r < DBL_MIN : giving up (-> +- Inf \n");
 - 
 
130 
#endif
 - 
 
131 
	    ML_ERROR(ME_RANGE);
 - 
 
132 
	    if(q < 0.0) return ML_NEGINF;
 - 
 
133 
	    else	return ML_POSINF;
 - 
 
134 
	}
 - 
 
135 
    }
 - 
 
136 
/* FIXME: This could be improved when log_p or !lower_tail ?
 - 
 
137 
 *	  (using p, not p_ , and a different derivative )
 - 
 
138 
 */
 - 
 
139 
#ifdef DEBUG_qnorm
 - 
 
140 
    REprintf("\t before final step: val = %7g\n", val);
 - 
 
141 
#endif
 - 
 
142 
    /* Final Newton step: */
 - 
 
143 
    val = val -
 - 
 
144 
	(pnorm(val, 0., 1., /*lower*/TRUE, /*log*/FALSE) - p_) /
 - 
 
145 
	 dnorm(val, 0., 1., /*log*/FALSE);
 - 
 
146 
 
 - 
 
147 
#else
 - 
 
148 
/*-- use AS 241 --- */
 76 
/*-- use AS 241 --- */
149 
/* double ppnd16_(double *p, long *ifault)*/
 77 
/* double ppnd16_(double *p, long *ifault)*/
150 
/*      ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
 78 
/*      ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
151 
 
 79 
 
152 
        Produces the normal deviate Z corresponding to a given lower
 80 
        Produces the normal deviate Z corresponding to a given lower
Line 217... Line 145...
217 
 
 145 
 
218 
	if(q < 0.0)
 146 
	if(q < 0.0)
219 
	    val = -val;
 147 
	    val = -val;
220 
        /* return (q >= 0.)? r : -r ;*/
 148 
        /* return (q >= 0.)? r : -r ;*/
221 
    }
 149 
    }
222 
 
 - 
 
223 
#endif
 - 
 
224 
/*-- Switch of AS 111 <-> AS 241 --- */
 - 
 
225 
 
 - 
 
226 
    return mu + sigma * val;
 150 
    return mu + sigma * val;
227 
}
 151 
}
228 
 
 152 
 
229 
 
 153 
 
230 
 
 154