The R Project SVN R

Rev

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

Rev 77685 Rev 88401
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
- 
 
 3 
 *  Copyright (C) 2000-2021 The R Core Team
3 
 *  Copyright (C) 1998 Ross Ihaka
 4 
 *  Copyright (C) 1998 Ross Ihaka
4 
 *  Copyright (C) 2000-12 The R Core Team
 - 
 
5 
 *
 5 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 6 
 *  This program is free software; you can redistribute it and/or modify
7 
 *  it under the terms of the GNU General Public License as published by
 7 
 *  it under the terms of the GNU General Public License as published by
8 
 *  the Free Software Foundation; either version 2 of the License, or
 8 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 9 
 *  (at your option) any later version.
Line 54... Line 54...
54 
 
 54 
 
55 
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)
56 
{
 56 
{
57 
    const static double eps = 1.e-15;
 57 
    const static double eps = 1.e-15;
58 
 
 58 
 
59 
    int kMax;
 - 
 
60 
    double k, ncp2, dx2, d, D;
 - 
 
61 
    LDOUBLE sum, term, p_k, q;
 - 
 
62 
 
 - 
 
63 
#ifdef IEEE_754
 59 
#ifdef IEEE_754
64 
    if (ISNAN(x) || ISNAN(a) || ISNAN(b) || ISNAN(ncp))
 60 
    if (ISNAN(x) || ISNAN(a) || ISNAN(b) || ISNAN(ncp))
65 
	return x + a + b + ncp;
 61 
	return x + a + b + ncp;
66 
#endif
 62 
#endif
67 
    if (ncp < 0 || a <= 0 || b <= 0)
 63 
    if (ncp < 0 || a <= 0 || b <= 0)
Line 72... Line 68...
72 
 
 68 
 
73 
    if (x < 0 || x > 1) return(R_D__0);
 69 
    if (x < 0 || x > 1) return(R_D__0);
74 
    if(ncp == 0)
 70 
    if(ncp == 0)
75 
	return dbeta(x, a, b, give_log);
 71 
	return dbeta(x, a, b, give_log);
76 
 
 72 
 
77 
    /* New algorithm, starting with *largest* term : */
 73 
    /* Non-central Beta:  New algorithm, starting with *largest* term : */
- 
 
 74 
    double
78 
    ncp2 = 0.5 * ncp;
 75 
	ncp2 = ldexp(ncp, -1), // = 0.5 * ncp
79 
    dx2 = ncp2*x;
 76 
	dx2 = ncp2*x,
80 
    d = (dx2 - a - 1)/2;
 77 
	d = ldexp(dx2 - a - 1, -1), // = (...)/2
81 
    D = d*d + dx2 * (a + b) - a;
 78 
	D = d*d + dx2 * (a + b) - a;
- 
 
 79 
    int kMax;
82 
    if(D <= 0) {
 80 
    if(D <= 0) {
83 
	kMax = 0;
 81 
	kMax = 0;
84 
    } else {
 82 
    } else {
85 
	D = ceil(d + sqrt(D));
 83 
	D = ceil(d + sqrt(D));
86 
	kMax = (D > 0) ? (int)D : 0;
 84 
	kMax = (D > 0) ? (int)D : 0;
87 
    }
 85 
    }
88 
 
 86 
 
- 
 
 87 
    LDOUBLE sum, term, p_k, q;
89 
    /* The starting "middle term" --- first look at it's log scale: */
 88 
    /* The starting "middle term" --- first look at it's log scale: */
90 
    term = dbeta(x, a + kMax, b, /* log = */ TRUE);
 89 
    term = dbeta(x, a + kMax, b, /* log = */ TRUE);
91 
    p_k = dpois_raw(kMax, ncp2,              TRUE);
 90 
    p_k = dpois_raw(kMax, ncp2,              TRUE);
92 
    if(x == 0. || !R_FINITE(term) || !R_FINITE((double)p_k)) /* if term = +Inf */
 91 
    if(x == 0. || !R_FINITE(term) || !R_FINITE((double)p_k)) /* if term = +Inf */
93 
	return R_D_exp((double)(p_k + term));
 92 
	return R_D_exp((double)(p_k + term));
Line 99... Line 98...
99 
    /* mid = 1 = the rescaled value, instead of  mid = exp(p_k); */
 98 
    /* mid = 1 = the rescaled value, instead of  mid = exp(p_k); */
100 
 
 99 
 
101 
    /* Now sum from the inside out */
 100 
    /* Now sum from the inside out */
102 
    sum = term = 1. /* = mid term */;
 101 
    sum = term = 1. /* = mid term */;
103 
    /* middle to the left */
 102 
    /* middle to the left */
104 
    k = kMax;
 103 
    double k = kMax;
105 
    while(k > 0 && term > sum * eps) {
 104 
    while(k > 0 && term > sum * eps) {
106 
	k--;
 105 
	k--;
107 
	q = /* 1 / r_k = */ (k+1)*(k+a) / (k+a+b) / dx2;
 106 
	q = /* 1 / r_k = */ (k+1)*(k+a) / (k+a+b) / dx2;
108 
	term *= q;
 107 
	term *= q;
109 
	sum += term;
 108 
	sum += term;