The R Project SVN R

Rev

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

Rev 44785 Rev 44932
Line 27... Line 27...
27 
 *    Computes the density of the noncentral beta distribution with
 27 
 *    Computes the density of the noncentral beta distribution with
28 
 *    noncentrality parameter ncp.  The noncentral beta distribution
 28 
 *    noncentrality parameter ncp.  The noncentral beta distribution
29 
 *    has density:
 29 
 *    has density:
30 
 *
 30 
 *
31 
 *		       Inf
 31 
 *		       Inf
32 
 *	  f(x|a,b,d) = SUM p(i) * B(a+i,b) * x^(a+i-1) * (1-x)^(b-1)
 32 
 *	f(x|a,b,ncp) = SUM  p(i) *  x^(a+i-1) * (1-x)^(b-1) / B(a+i,b)
33 
 *		       i=0
 33 
 *		       i=0
34 
 *
 34 
 *
35 
 *    where:
 35 
 *    where:
36 
 *
 36 
 *
37 
 *		p(k) = exp(-ncp) ncp^k / k!
 37 
 *		p(k) = exp(-ncp/2) (ncp/2)^k / k!
38 
 *
 38 
 *
39 
 *	      B(a,b) = Gamma(a+b) / (Gamma(a) * Gamma(b))
 39 
 *	      B(a,b) = Gamma(a) * Gamma(b) / Gamma(a+b)
40 
 *
 40 
 *
41 
 *
 41 
 *
42 
 *    This can be computed efficiently by using the recursions:
 42 
 *    This can be computed efficiently by using the recursions:
43 
 *
 43 
 *
44 
 *	      p(k+1) = (ncp/(k+1)) * p(k-1)
 44 
 *	      p(k+1) = ncp/2 / (k+1) * p(k)
45 
 *
 45 
 *
46 
 *	  B(a+k+1,b) = ((a+b+k)/(a+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
 48 
 *    The summation of the series continues until
49 
 *
 49 
 *
50 
 *		psum = p(0) + ... + p(k)
 50 
 *		psum = p(0) + ... + p(k)
51 
 *
 51 
 *
Line 57... Line 57...
57 
#include "dpq.h"
 57 
#include "dpq.h"
58 
 
 58 
 
59 
double dnbeta(double x, double a, double b, double ncp, int give_log)
 59 
double dnbeta(double x, double a, double b, double ncp, int give_log)
60 
{
 60 
{
61 
    const static double eps = 1.e-14;
 61 
    const static double eps = 1.e-14;
62 
    const int maxiter = 200;
 62 
    const int maxiter = 10000; /* was 200 */
63 
 
 63 
 
64 
    double k, ncp2;
 64 
    double k, ncp2;
65 
    LDOUBLE psum, sum, term, weight;
 65 
    LDOUBLE psum, sum, term, weight;
66 
 
 66 
 
67 
#ifdef IEEE_754
 67 
#ifdef IEEE_754