The R Project SVN R

Rev

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

Rev 44946 Rev 45272
Line 43... Line 43...
43 
 *
 43 
 *
44 
 *	      p(k+1) = ncp/2 / (k+1) * p(k)
 44 
 *	      p(k+1) = ncp/2 / (k+1) * p(k)
45 
 *
 45 
 *
46 
 *      B(a+k+1,b)   = (a+k)/(a+b+k) * B(a+k,b)
 46 
 *      B(a+k+1,b)   = (a+k)/(a+b+k) * B(a+k,b)
47 
 *
 47 
 *
48 
 *    The summation of the series continues until
 - 
 
49 
 *
 - 
 
50 
 *		psum = p(0) + ... + p(k)
 - 
 
51 
 *
 - 
 
52 
 *    is close to 1.  Here we continue until 1 - psum < epsilon,
 48 
 * The new algorithm first determines for which k the k-th term is maximal,
53 
 *    with epsilon set close to the relative machine precision.
 49 
 * and then sums outwards to both sides from the 'mid'.
54 
 */
 50 
 */
55 
 
 51 
 
56 
#include "nmath.h"
 52 
#include "nmath.h"
57 
#include "dpq.h"
 53 
#include "dpq.h"
58 
 
 54 
 
59 
double dnbeta(double x, double a, double b, double ncp, int give_log)
 55 
double dnbeta(double x, double a, double b, double ncp, int give_log)
60 
{
 56 
{
61 
    const static double eps = 1.e-14;
 57 
    const static double eps = 1.e-15;
62 
    const int maxiter = 10000; /* was 200 */
 - 
 
63 
 
 58 
 
- 
 
 59 
    int kMax;
64 
    double k, ncp2;
 60 
    double k, ncp2, dx2, d, D;
65 
    LDOUBLE psum, sum, term, weight;
 61 
    LDOUBLE sum, term, p_k, q;
66 
 
 62 
 
67 
#ifdef IEEE_754
 63 
#ifdef IEEE_754
68 
    if (ISNAN(x) || ISNAN(a) || ISNAN(b) || ISNAN(ncp))
 64 
    if (ISNAN(x) || ISNAN(a) || ISNAN(b) || ISNAN(ncp))
69 
	return x + a + b + ncp;
 65 
	return x + a + b + ncp;
70 
#endif
 66 
#endif
Line 76... Line 72...
76 
 
 72 
 
77 
    if (x < 0 || x > 1) return(R_D__0);
 73 
    if (x < 0 || x > 1) return(R_D__0);
78 
    if(ncp == 0)
 74 
    if(ncp == 0)
79 
	return dbeta(x, a, b, give_log);
 75 
	return dbeta(x, a, b, give_log);
80 
 
 76 
 
81 
    term = dbeta(x, a, b, /* log = */ TRUE);
 - 
 
82 
    if(!R_FINITE(term)) /* in particular, if term = +Inf */
 77 
    /* New algorithm, starting with *largest* term : */
83 
	return R_D_exp(term);
 - 
 
84 
    ncp2 = 0.5 * ncp;
 78 
    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);
 79 
    dx2 = ncp2*x;
88 
    weight = exp(- ncp2);
 80 
    d = (dx2 - a - 1)/2;
89 
    sum = weight * term;
 81 
    D = d*d + dx2 * (a + b) - a;
90 
    if(sum == 0.) {
 82 
    if(D <= 0) {
91 
	if(term != 0.) /* (x = {0,1} gives true 0 for a,b>=1) */
 - 
 
92 
	    ML_ERROR(ME_UNDERFLOW, "dnbeta");
 83 
	kMax = 0;
93 
    } else {
 84 
    } else {
94 
	psum = weight;
 - 
 
95 
	for(k = 1; k <= maxiter; k++) {
 - 
 
96 
	    double c1, c2, t;
 - 
 
97 
	    weight *= (c1 = (ncp2 / k));
 85 
	D = ceil(d + sqrt(D));
98 
	    term   *= (c2 = x * (a + b) / a);
 - 
 
99 
	    sum  += (t = weight * term);
 86 
	kMax = (D > 0) ? (int)D : 0;
100 
	    psum += weight;
 - 
 
101 
	    a += 1;
 - 
 
102 
	    if(c1*c2 < 1 && psum + eps > 1 && t < eps * sum)
 - 
 
103 
		break;
 - 
 
104 
	    else if(k == maxiter) /* not converged */
 - 
 
105 
		ML_ERROR(ME_NOCONV, "dnbeta");
 - 
 
106 
	}
 - 
 
107 
    }
 87 
    }
- 
 
 88 
 
- 
 
 89 
    /* The starting "middle term" --- first look at it's log scale: */
- 
 
 90 
    term = dbeta(x, a + kMax, b, /* log = */ TRUE);
- 
 
 91 
    p_k = dpois_raw(kMax, ncp2,              TRUE);
- 
 
 92 
    if(x == 0. || !R_FINITE(term) || !R_FINITE(p_k)) /* if term = +Inf */
- 
 
 93 
	return R_D_exp(p_k + term);
- 
 
 94 
 
- 
 
 95 
    /* Now if s_k := p_k * t_k  {here = exp(p_k + term)} would underflow,
- 
 
 96 
     * we should rather scale everything and re-scale at the end:*/
- 
 
 97 
 
- 
 
 98 
    p_k += term; /* = log(p_k) + log(t_k) == log(s_k) -- used at end to rescale */
- 
 
 99 
    /* mid = 1 = the rescaled value, instead of  mid = exp(p_k); */
- 
 
 100 
 
- 
 
 101 
    /* Now sum from the inside out */
- 
 
 102 
    sum = term = 1. /* = mid term */;
- 
 
 103 
    /* middle to the left */
- 
 
 104 
    k = kMax;
- 
 
 105 
    while(k > 0 && term > sum * eps) {
- 
 
 106 
	k--;
- 
 
 107 
	q = /* 1 / r_k = */ (k+1)*(k+a) / (k+a+b) / dx2;
- 
 
 108 
	term *= q;
- 
 
 109 
	sum += term;
- 
 
 110 
    }
- 
 
 111 
    /* middle to the right */
- 
 
 112 
    term = 1.;
- 
 
 113 
    k = kMax;
- 
 
 114 
    do {
- 
 
 115 
	q = /* r_{old k} = */ dx2 * (k+a+b) / (k+a) / (k+1);
- 
 
 116 
	k++;
- 
 
 117 
	term *= q;
- 
 
 118 
	sum += term;
- 
 
 119 
    } while (term > sum * eps);
- 
 
 120 
 
108 
    return R_D_val(sum);
 121 
    return R_D_exp(p_k + log(sum));
109 
}
 122 
}