The R Project SVN R

Rev

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

Rev 44932 Rev 44946
Line 76... Line 76...
76 
 
 76 
 
77 
    if (x < 0 || x > 1) return(R_D__0);
 77 
    if (x < 0 || x > 1) return(R_D__0);
78 
    if(ncp == 0)
 78 
    if(ncp == 0)
79 
	return dbeta(x, a, b, give_log);
 79 
	return dbeta(x, a, b, give_log);
80 
 
 80 
 
81 
    term = dbeta(x, a, b, /* log = */ FALSE);
 81 
    term = dbeta(x, a, b, /* log = */ TRUE);
82 
    if(!R_FINITE(term)) /* in particular, if term = +Inf */
 82 
    if(!R_FINITE(term)) /* in particular, if term = +Inf */
83 
	return R_D_val(term);
 83 
	return R_D_exp(term);
84 
    ncp2 = 0.5 * ncp;
 84 
    ncp2 = 0.5 * ncp;
- 
 
 85 
    /* FIXME: prevent underflow in term *and* weight-- probably should
- 
 
 86 
     * work in log scale and *rescale* when needed ..*/
- 
 
 87 
    term   = exp(term);
85 
    weight = exp(- ncp2);
 88 
    weight = exp(- ncp2);
86 
    sum	 = weight * term;
 89 
    sum = weight * term;
- 
 
 90 
    if(sum == 0.) {
- 
 
 91 
	if(term != 0.) /* (x = {0,1} gives true 0 for a,b>=1) */
- 
 
 92 
	    ML_ERROR(ME_UNDERFLOW, "dnbeta");
- 
 
 93 
    } else {
87 
    psum = weight;
 94 
	psum = weight;
88 
    for(k = 1; k <= maxiter; k++) {
 95 
	for(k = 1; k <= maxiter; k++) {
- 
 
 96 
	    double c1, c2, t;
89 
	weight *= (ncp2 / k);
 97 
	    weight *= (c1 = (ncp2 / k));
90 
	term *= x * (a + b) / a;
 98 
	    term   *= (c2 = x * (a + b) / a);
91 
	sum  += weight * term;
 99 
	    sum  += (t = weight * term);
92 
	psum += weight;
 100 
	    psum += weight;
93 
	a += 1;
 101 
	    a += 1;
94 
	if(1 - psum < eps) break;
 102 
	    if(c1*c2 < 1 && psum + eps > 1 && t < eps * sum)
95 
    }
 103 
		break;
96 
    if(1 - psum >= eps) { /* not converged */
 104 
	    else if(k == maxiter) /* not converged */
97 
	ML_ERROR(ME_PRECISION, "dnbeta");
 105 
		ML_ERROR(ME_NOCONV, "dnbeta");
- 
 
 106 
	}
98 
    }
 107 
    }
99 
    return R_D_val(sum);
 108 
    return R_D_val(sum);
100 
}
 109 
}