The R Project SVN R

Rev

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

Rev 87903 Rev 88401
Line 11... Line 11...
11 
 *  Returns the cumulative probability of x for the non-central
 11 
 *  Returns the cumulative probability of x for the non-central
12 
 *  beta distribution with parameters a, b and non-centrality ncp.
 12 
 *  beta distribution with parameters a, b and non-centrality ncp.
13 
 *
 13 
 *
14 
 *  Auxiliary routines required:
 14 
 *  Auxiliary routines required:
15 
 *	lgamma - log-gamma function
 15 
 *	lgamma - log-gamma function
16 
 *      pbeta  - incomplete-beta function {nowadays: pbeta_raw() -> bratio()}
 16 
 *      pbeta  - incomplete-beta function {nowadays: directly bratio()}
17 
 */
 17 
 */
18 
 
 18 
 
19 
#include "nmath.h"
 19 
#include "nmath.h"
20 
#include "dpq.h"
 20 
#include "dpq.h"
21 
 
 21 
 
Line 47... Line 47...
47 
    x0 = floor(fmax2(c - 7. * sqrt(c), 0.));
 47 
    x0 = floor(fmax2(c - 7. * sqrt(c), 0.));
48 
    a0 = a + x0;
 48 
    a0 = a + x0;
49 
    lBeta = lbeta(a0, b); // = lgammafn(a0) + lgammafn(b) - lgammafn(a0 + b);
 49 
    lBeta = lbeta(a0, b); // = lgammafn(a0) + lgammafn(b) - lgammafn(a0 + b);
50 
    /* temp = pbeta_raw(x, a0, b, TRUE, FALSE), but using (x, o_x): */
 50 
    /* temp = pbeta_raw(x, a0, b, TRUE, FALSE), but using (x, o_x): */
51 
    bratio(a0, b, x, o_x, &temp, &tmp_c, &ierr, FALSE);
 51 
    bratio(a0, b, x, o_x, &temp, &tmp_c, &ierr, FALSE);
- 
 
 52 
    // TODO(?): report error code as in  ./pbeta.c
52 
 
 53 
 
53 
    gx = exp(a0 * log(x) + b * (x < .5 ? log1p(-x) : log(o_x))
 54 
    gx = exp(a0 * log(x) + b * (x < .5 ? log1p(-x) : log(o_x))
54 
	     - lBeta - log(a0));
 55 
	     - lBeta - log(a0));
55 
    if (a0 > a) // x0 >= 1 (and *not* x0 << a)
 56 
    if (a0 > a) // x0 >= 1 (and *not* x0 << a)
56 
	q = exp(-c + x0 * log(c) - lgammafn(x0 + 1.));
 57 
	q = exp(-c + x0 * log(c) - lgammafn(x0 + 1.));