The R Project SVN R

Rev

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

Rev 42302 Rev 42303
Line -... Line 1...
- 
 
 1 
/*
- 
 
 2 
 *  AUTHOR
- 
 
 3 
 *    Claus Ekstrøm, ekstrom@dina.kvl.dk
- 
 
 4 
 *    July 15, 2003.
- 
 
 5 
 *
- 
 
 6 
 *  Merge in to R:
- 
 
 7 
 *	Copyright (C) 2003 The R Foundation
- 
 
 8 
 *
- 
 
 9 
 *  This program is free software; you can redistribute it and/or modify
- 
 
 10 
 *  it under the terms of the GNU General Public License as published by
- 
 
 11 
 *  the Free Software Foundation; either version 2 of the License, or
- 
 
 12 
 *  (at your option) any later version.
- 
 
 13 
 *
- 
 
 14 
 *  This program is distributed in the hope that it will be useful,
- 
 
 15 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- 
 
 16 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- 
 
 17 
 *  GNU General Public License for more details.
- 
 
 18 
 *
- 
 
 19 
 *  You should have received a copy of the GNU General Public License
- 
 
 20 
 *  along with this program; if not, a copy is available at
- 
 
 21 
 *  http://www.r-project.org/Licenses/
- 
 
 22 
 *
- 
 
 23 
 *
- 
 
 24 
 *  NOTE
- 
 
 25 
 *
- 
 
 26 
 *    Requires the following auxiliary routines:
- 
 
 27 
 *
- 
 
 28 
 *	lgammafn(x)	- log gamma function
- 
 
 29 
 *	pnt(x, df, ncp) - the distribution function for
- 
 
 30 
 *			  the non-central t distribution
- 
 
 31 
 *
- 
 
 32 
 *
- 
 
 33 
 * DESCRIPTION
- 
 
 34 
 *
- 
 
 35 
 *    The non-central t density is
- 
 
 36 
 *
- 
 
 37 
 *	   f(x, df, ncp) =
- 
 
 38 
 *		df^(df/2) * exp(-.5*ncp^2) /
- 
 
 39 
 *		(sqrt(pi)*gamma(df/2)*(df+x^2)^((df+1)/2)) *
- 
 
 40 
 *		sum_{k=0}^Inf  gamma((df + k + df)/2)*ncp^k /
- 
 
 41 
 *				prod(1:k)*(2*x^2/(df+x^2))^(k/2)
- 
 
 42 
 *
- 
 
 43 
 *    The functional relationship
- 
 
 44 
 *
- 
 
 45 
 *	   f(x, df, ncp) = df/x *
- 
 
 46 
 *			      (F(sqrt((df+2)/df)*x, df+2, ncp) - F(x, df, ncp))
- 
 
 47 
 *
- 
 
 48 
 *    is used to evaluate the density at x != 0 and
- 
 
 49 
 *
- 
 
 50 
 *	   f(0, df, ncp) = exp(-.5*ncp^2) /
- 
 
 51 
 *				(sqrt(pi)*sqrt(df)*gamma(df/2))*gamma((df+1)/2)
- 
 
 52 
 *
- 
 
 53 
 *    is used for x=0.
- 
 
 54 
 *
- 
 
 55 
 *    All calculations are done on log-scale to increase stability.
- 
 
 56 
 *
- 
 
 57 
 */
- 
 
 58 
 
- 
 
 59 
#include "nmath.h"
- 
 
 60 
#include "dpq.h"
- 
 
 61 
 
- 
 
 62 
double dnt(double x, double df, double ncp, int give_log)
- 
 
 63 
{
- 
 
 64 
    double u;
- 
 
 65 
#ifdef IEEE_754
- 
 
 66 
    if (ISNAN(x) || ISNAN(df))
- 
 
 67 
	return x + df;
- 
 
 68 
#endif
- 
 
 69 
 
- 
 
 70 
    /* If non-positive df then error */
- 
 
 71 
    if (df <= 0.0) ML_ERR_return_NAN;
- 
 
 72 
 
- 
 
 73 
    if(ncp == 0.0) return dt(x, df, give_log);
- 
 
 74 
 
- 
 
 75 
    /* If x is infinite then return 0 */
- 
 
 76 
    if(!R_FINITE(x))
- 
 
 77 
	return R_D__0;
- 
 
 78 
 
- 
 
 79 
    /* If infinite df then the density is identical to a
- 
 
 80 
       normal distribution with mean = ncp.  However, the formula
- 
 
 81 
       loses a lot of accuracy around df=1e9
- 
 
 82 
    */
- 
 
 83 
    if(!R_FINITE(df) || df > 1e8)
- 
 
 84 
	return dnorm(x, ncp, 1., give_log);
- 
 
 85 
 
- 
 
 86 
    /* Do calculations on log scale to stabilize */
- 
 
 87 
 
- 
 
 88 
    /* Consider two cases: x ~= 0 or not */
- 
 
 89 
    if (fabs(x) > sqrt(df * DBL_EPSILON)) {
- 
 
 90 
	u = log(df) - log(fabs(x)) +
- 
 
 91 
	    log(fabs(pnt(x*sqrt((df+2)/df), df+2, ncp, 1, 0) -
- 
 
 92 
		     pnt(x, df, ncp, 1, 0)));
- 
 
 93 
	/* FIXME: the above still suffers from cancellation (but not horribly) */
- 
 
 94 
    }
- 
 
 95 
    else {  /* x ~= 0 : -> same value as for  x = 0 */
- 
 
 96 
	u = lgammafn((df+1)/2) - lgammafn(df/2)
- 
 
 97 
	    - .5*(log(M_PI) + log(df) + ncp*ncp);
- 
 
 98 
    }
- 
 
 99 
 
- 
 
 100 
    return (give_log ? u : exp(u));
- 
 
 101 
}