The R Project SVN R

Rev

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

Rev 59037 Rev 61551
Line 35... Line 35...
35 
#include "nmath.h"
 35 
#include "nmath.h"
36 
#include "dpq.h"
 36 
#include "dpq.h"
37 
 
 37 
 
38 
double dt(double x, double n, int give_log)
 38 
double dt(double x, double n, int give_log)
39 
{
 39 
{
40 
    double t, u;
 - 
 
41 
#ifdef IEEE_754
 40 
#ifdef IEEE_754
42 
    if (ISNAN(x) || ISNAN(n))
 41 
    if (ISNAN(x) || ISNAN(n))
43 
	return x + n;
 42 
	return x + n;
44 
#endif
 43 
#endif
45 
 
 - 
 
46 
    if (n <= 0) ML_ERR_return_NAN;
 44 
    if (n <= 0) ML_ERR_return_NAN;
47 
    if(!R_FINITE(x))
 45 
    if(!R_FINITE(x))
48 
	return R_D__0;
 46 
	return R_D__0;
49 
    if(!R_FINITE(n))
 47 
    if(!R_FINITE(n))
50 
	return dnorm(x, 0., 1., give_log);
 48 
	return dnorm(x, 0., 1., give_log);
51 
 
 49 
 
52 
    t = -bd0(n/2.,(n+1)/2.) + stirlerr((n+1)/2.) - stirlerr(n/2.);
 50 
    double u, ax, t = -bd0(n/2.,(n+1)/2.) + stirlerr((n+1)/2.) - stirlerr(n/2.),
- 
 
 51 
	x2n = x*x/n, // in  [0, Inf]
- 
 
 52 
	l_x2n; // := log(sqrt(1 + x2n)) = log(1 + x2n)/2
- 
 
 53 
    Rboolean lrg_x2n =  (x2n > 1./DBL_EPSILON);
- 
 
 54 
    if (lrg_x2n) { // large x^2/n :
- 
 
 55 
	ax = fabs(x);
- 
 
 56 
	l_x2n = log(ax) - log(n)/2.; // = log(x2n)/2 = 1/2 * log(x^2 / n)
- 
 
 57 
	u = //  log(1 + x2n) * n/2 =  n * log(1 + x2n)/2 =
- 
 
 58 
	    n * l_x2n;
- 
 
 59 
    }
53 
    if ( x*x > 0.2*n )
 60 
    else if (x2n > 0.2) {
54 
	u = log( 1+ x*x/n ) * n/2;
 61 
	l_x2n = log(1 + x2n)/2.;
- 
 
 62 
	u = n * l_x2n;
55 
    else
 63 
    } else {
- 
 
 64 
	l_x2n = log1p(x2n)/2.;
56 
	u = -bd0(n/2.,(n+x*x)/2.) + x*x/2.;
 65 
	u = -bd0(n/2.,(n+x*x)/2.) + x*x/2.;
- 
 
 66 
    }
- 
 
 67 
 
- 
 
 68 
    //old: return  R_D_fexp(M_2PI*(1+x2n), t-u);
- 
 
 69 
 
- 
 
 70 
    // R_D_fexp(f,x) :=  (give_log ? -0.5*log(f)+(x) : exp(x)/sqrt(f))
- 
 
 71 
    // f = 2pi*(1+x2n)
- 
 
 72 
    //  ==> 0.5*log(f) = log(2pi)/2 + log(1+x2n)/2 = log(2pi)/2 + l_x2n
- 
 
 73 
    //	     1/sqrt(f) = 1/sqrt(2pi * (1+ x^2 / n))
- 
 
 74 
    //		       = 1/sqrt(2pi)/(|x|/sqrt(n)*sqrt(1+1/x2n))
- 
 
 75 
    //		       = M_1_SQRT_2PI * sqrt(n)/ (|x|*sqrt(1+1/x2n))
- 
 
 76 
    if(give_log)
- 
 
 77 
	return t-u - (M_LN_SQRT_2PI + l_x2n);
57 
 
 78 
 
- 
 
 79 
    // else :  if(lrg_x2n) : sqrt(1 + 1/x2n) ='= sqrt(1) = 1
- 
 
 80 
    double I_sqrt_ = (lrg_x2n ? sqrt(n)/ax : exp(-l_x2n));
58 
    return R_D_fexp(M_2PI*(1+x*x/n), t-u);
 81 
    return exp(t-u) * M_1_SQRT_2PI * I_sqrt_;
59 
}
 82 
}