The R Project SVN R

Rev

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

Rev 11723 Rev 14563
Line 38... Line 38...
38 
#include "dpq.h"
 38 
#include "dpq.h"
39 
 
 39 
 
40 
double dbeta(double x, double a, double b, int give_log)
 40 
double dbeta(double x, double a, double b, int give_log)
41 
{
41 
{
42 
    double f, p;
 42 
    double f, p;
- 
 
 43 
    volatile double am1, bm1; /* prevent roundoff trouble on some
- 
 
 44 
                                 platforms */
43 
 
 45 
 
44 
#ifdef IEEE_754
 46 
#ifdef IEEE_754
45 
    /* NaNs propagated correctly */
 47 
    /* NaNs propagated correctly */
46 
    if (ISNAN(x) || ISNAN(a) || ISNAN(b)) return x + a + b;
 48 
    if (ISNAN(x) || ISNAN(a) || ISNAN(b)) return x + a + b;
47 
#endif
 49 
#endif
Line 63... Line 65...
63 
	    f = a*b/((a+b)*x*(1-x));
 65 
	    f = a*b/((a+b)*x*(1-x));
64 
	    p = dbinom_raw(a,a+b, x,1-x, give_log);
 66 
	    p = dbinom_raw(a,a+b, x,1-x, give_log);
65 
	}
 67 
	}
66 
	else {			/* a < 1 <= b */
 68 
	else {			/* a < 1 <= b */
67 
	    f = a/x;
 69 
	    f = a/x;
- 
 
 70 
	    bm1 = b - 1;
68 
	    p = dbinom_raw(a,a+b-1, x,1-x, give_log);
 71 
	    p = dbinom_raw(a,a+bm1, x,1-x, give_log);
69 
	}
 72 
	}
70 
    }
 73 
    }
71 
    else {
74 
    else {
72 
	if (b < 1) {		/* a >= 1 > b */
 75 
	if (b < 1) {		/* a >= 1 > b */
73 
	    f = b/(1-x);
 76 
	    f = b/(1-x);
- 
 
 77 
	    am1 = a - 1;
74 
	    p = dbinom_raw(a-1,a+b-1, x,1-x, give_log);
 78 
	    p = dbinom_raw(am1,am1+b, x,1-x, give_log);
75 
	}
 79 
	}
76 
	else {			/* a,b >= 1 */
 80 
	else {			/* a,b >= 1 */
77 
	    f = a+b-1;
 81 
	    f = a+b-1;
- 
 
 82 
	    am1 = a - 1;
- 
 
 83 
	    bm1 = b - 1;
78 
	    p = dbinom_raw(a-1,a+b-2, x,1-x, give_log);
 84 
	    p = dbinom_raw(am1,am1+bm1, x,1-x, give_log);
79 
	}
 85 
	}
80 
    }
 86 
    }
81 
    return( (give_log) ? p + log(f) : p*f );
 87 
    return( (give_log) ? p + log(f) : p*f );
82 
}
 88 
}