The R Project SVN R

Rev

Rev 82839 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 82839 Rev 88401
Line 79... Line 79...
79 
    if(!R_FINITE(x))
 79 
    if(!R_FINITE(x))
80 
	return R_D__0;
 80 
	return R_D__0;
81 
 
 81 
 
82 
    /* If infinite df then the density is identical to a
 82 
    /* If infinite df then the density is identical to a
83 
       normal distribution with mean = ncp.  However, the formula
 83 
       normal distribution with mean = ncp.  However, the formula
84 
       loses a lot of accuracy around df=1e9
 84 
       loses a lot of accuracy around df=1e9 // FIXME?
85 
    */
 85 
    */
86 
    if(!R_FINITE(df) || df > 1e8)
 86 
    if(!R_FINITE(df) || df > 1e8)
87 
	return dnorm(x, ncp, 1., give_log);
 87 
	return dnorm(x, ncp, 1., give_log);
88 
 
 88 
 
89 
    /* Do calculations on log scale to stabilize */
 89 
    /* Do calculations on log scale to stabilize */
90 
 
 90 
 
91 
    /* Consider two cases: x ~= 0 or not */
 91 
    /* Consider two cases: x ~= 0 or not */
92 
    if (fabs(x) > sqrt(df * DBL_EPSILON)) {
 92 
    if (fabs(x) > sqrt(df * DBL_EPSILON)) { // |x| > eps * sqrt(df)
93 
	u = log(df) - log(fabs(x)) +
 93 
	u = log(df) - log(fabs(x)) +
94 
	    log(fabs(pnt(x*sqrt((df+2)/df), df+2, ncp, 1, 0) -
 94 
	    log(fabs(pnt(x*sqrt((df+2)/df), df+2, ncp, 1, 0) -
95 
		     pnt(x, df, ncp, 1, 0)));
 95 
		     pnt(x, df, ncp, 1, 0)));
96 
	/* FIXME: the above still suffers from cancellation (but not horribly) */
 96 
	/* FIXME: the above still suffers from cancellation (but not horribly) */
97 
    }
 97 
    }