The R Project SVN R

Rev

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

Rev 13115 Rev 24049
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Algorithm AS 275 Appl.Statist. (1992), vol.41, no.2
 2 
 *  Algorithm AS 275 Appl.Statist. (1992), vol.41, no.2
3 
 *  original  (C) 1992 Royal Statistical Society
 3 
 *  original  (C) 1992	     Royal Statistical Society
4 
 *  Copyright (C) 2000 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 
 *  Computes the noncentral chi-squared distribution function with
 7 
 *  Computes the noncentral chi-squared distribution function with
7 
 *  positive real degrees of freedom f and nonnegative noncentrality
 8 
 *  positive real degrees of freedom f and nonnegative noncentrality
8 
 *  parameter theta
 9 
 *  parameter theta
9 
 */
 10 
 */
Line 31... Line 32...
31 
    if (f < 0. || theta < 0.) ML_ERR_return_NAN;
 32 
    if (f < 0. || theta < 0.) ML_ERR_return_NAN;
32 
 
 33 
 
33 
    return (R_DT_val(pnchisq_raw(x, f, theta, 1e-12, 10000)));
 34 
    return (R_DT_val(pnchisq_raw(x, f, theta, 1e-12, 10000)));
34 
}
 35 
}
35 
 
 36 
 
36 
double pnchisq_raw(double x, double f, double theta,
37 
double pnchisq_raw(double x, double f, double theta,
37 
		   double errmax, int itrmax)
 38 
		   double errmax, int itrmax)
38 
{
 39 
{
39 
    double ans, lam, u, v, x2, f2, t, term, bound, f_x_2n, f_2n;
 40 
    double ans, lam, u, v, x2, f2, t, term, bound, f_x_2n, f_2n;
40 
    int n, flag;
 41 
    int n, flag;
41 
 
 42 
 
- 
 
 43 
    static const double my_dbl_min_exp = M_LN2 * DBL_MIN_EXP;
- 
 
 44 
    /*= -708.3964 for IEEE double precision */
- 
 
 45 
 
42 
    if (x <= 0.)	return 0.;
 46 
    if (x <= 0.)	return 0.;
43 
    if(!R_FINITE(x))	return 1.;
 47 
    if(!R_FINITE(x))	return 1.;
44 
 
 48 
 
45 
    lam = .5 * theta;
 - 
 
46 
#ifdef DEBUG_pnch
 49 
#ifdef DEBUG_pnch
47 
    REprintf("pnchisq(x=%12g, f=%12g, theta=%12g):\n",x,f,theta);
 50 
    REprintf("pnchisq(x=%g, f=%g, theta=%g): ",x,f,theta);
48 
#endif
 51 
#endif
- 
 
 52 
    lam = .5 * theta;
- 
 
 53 
    if(-lam < my_dbl_min_exp)
- 
 
 54 
	MATHLIB_ERROR("non centrality parameter (= %g) too large for current algorithm", theta)
49 
 
 55 
 
50 
    /* evaluate the first term */
 56 
    /* evaluate the first term */
51 
 
 57 
 
52 
    v = u = exp(-lam);
 58 
    v = u = exp(-lam);
53 
    x2 = .5 * x;
 59 
    x2 = .5 * x;
54 
    f2 = .5 * f;
 60 
    f2 = .5 * f;
55 
    f_x_2n = f - x;
 61 
    f_x_2n = f - x;
- 
 
 62 
 
- 
 
 63 
#ifdef DEBUG_pnch
- 
 
 64 
    REprintf("-- v=exp(-th/2)=%g, x/2= %g, f/2= %g\n",v,x2,f2);
- 
 
 65 
#endif
- 
 
 66 
 
56 
    if(f2 * DBL_EPSILON > 0.125 &&
 67 
    if(f2 * DBL_EPSILON > 0.125 &&
57 
       fabs(t = x2 - f2) < sqrt(DBL_EPSILON) * f2) {
 68 
       fabs(t = x2 - f2) < sqrt(DBL_EPSILON) * f2) {
58 
	/* evade cancellation error */
 69 
	/* evade cancellation error */
59 
	t = exp((1 - t)*(2 - t/(f2 + 1))) / sqrt(2*M_PI*(f2 + 1));
 70 
	t = exp((1 - t)*(2 - t/(f2 + 1))) / sqrt(2*M_PI*(f2 + 1));
- 
 
 71 
#ifdef DEBUG_pnch
- 
 
 72 
	REprintf(" (case I) ==> t= %g\n",t);
- 
 
 73 
#endif
60 
    }
 74 
    }
61 
    else {
 75 
    else {
62 
	/* careful not to overflow .. : */
 76 
	/* careful not to overflow .. : */
63 
	t = exp(f2*log(x2) -x2 - lgammafn(f2 + 1));
 77 
	t = f2*log(x2) -x2 - lgammafn(f2 + 1);
- 
 
 78 
	if (t < my_dbl_min_exp &&
- 
 
 79 
	    x > f + theta +  3* sqrt( 2*(f + 2*theta))) {
- 
 
 80 
	    /* x > E[X] + 3* sigma(X) */
- 
 
 81 
	    return 1.; /* better than 0 ! */
- 
 
 82 
	} /* else */
- 
 
 83 
	t = exp(t);
- 
 
 84 
#ifdef DEBUG_pnch
- 
 
 85 
	REprintf(" (case 2) ==> t= %g\n",t);
- 
 
 86 
#endif
64 
    }
 87 
    }
65 
 
 88 
 
- 
 
 89 
    if(t <= 0)
- 
 
 90 
	MATHLIB_WARNING3("too large x (= %g) or centrality parameter %g for current algorithm;\n\t%s", theta, x, "result is probably invalid!");
- 
 
 91 
 
66 
    term = v * t;
 92 
    term = v * t;
67 
    ans = term;
 93 
    ans = term;
68 
#ifdef DEBUG_pnch
 - 
 
69 
    REprintf("\t v=exp(-th/2)=%12g, x/2=%12g, f/2=%12g ==> t=%12g\n",v,x2,f2,t);
 - 
 
70 
#endif
 - 
 
71 
 
 94 
 
72 
    /* check if (f+2n) is greater than x */
 95 
    /* check if (f+2n) is greater than x */
73 
 
 96 
 
74 
    flag = FALSE;
 97 
    flag = FALSE;
75 
    n = 1;
 98 
    n = 1;
76 
    f_2n = f + 2.;/* = f + 2*n */
 99 
    f_2n = f + 2.;/* = f + 2*n */
77 
    f_x_2n += 2.;/* = f - x + 2*n */
 100 
    f_x_2n += 2.;/* = f - x + 2*n */
78 
    for(;;) {
 101 
    for(;;) {
79 
#ifdef DEBUG_pnch
 102 
#ifdef DEBUG_pnch
80 
	REprintf(" _OL_: n=%d",n);
 103 
	REprintf("\n _OL_: n=%d",n);
81 
#endif
 104 
#endif
82 
	if (f_x_2n > 0) {
 105 
	if (f_x_2n > 0) {
83 
 
 106 
 
84 
	    /* find the error bound and check for convergence */
 107 
	    /* find the error bound and check for convergence */
85 
	    flag = TRUE;
 108 
	    flag = TRUE;
Line 101... Line 124...
101 
	    if (!flag && n <= itrmax)
 124 
	    if (!flag && n <= itrmax)
102 
		break;
 125 
		break;
103 
	L10:
 126 
	L10:
104 
	    bound = t * x / f_x_2n;
 127 
	    bound = t * x / f_x_2n;
105 
#ifdef DEBUG_pnch
 128 
#ifdef DEBUG_pnch
106 
	    REprintf("\tL10: n=%d; term=%12g; bound=%12g\n",n,term,bound);
 129 
	    REprintf("\n L10: n=%d; term= %g; bound= %g",n,term,bound);
107 
#endif
 130 
#endif
108 
	    if (bound <= errmax || n > itrmax)
 131 
	    if (bound <= errmax || n > itrmax)
109 
		goto L_End;
 132 
		goto L_End;
110 
	}
 133 
	}
111 
    }
 134 
    }
112 
L_End:
 135 
L_End:
113 
    if (bound > errmax) { /* NOT converged */
 136 
    if (bound > errmax) { /* NOT converged */
114 
	ML_ERROR(ME_PRECISION);
 137 
	ML_ERROR(ME_PRECISION);
115 
    }
 138 
    }
116 
#ifdef DEBUG_pnch
 139 
#ifdef DEBUG_pnch
117 
    REprintf("\tL_End: n=%d; term=%12g; bound=%12g\n",n,term,bound);
 140 
    REprintf("\n == L_End: n=%d; term= %g; bound=%g\n",n,term,bound);
118 
#endif
 141 
#endif
119 
    return (ans);
 142 
    return (ans);
120 
}
 143 
}