[R] Density of non-central t distribution

[Wolfgang Viechtbauer]

For those who might be interested:

I wrote a function that gives the pdf of a non-central t distribution.
This is an approximation based on Resnikoff & Lieberman (1957):

m   = degrees of freedom
ncp = non-centrality parameter

dtnc <- function(x, m, ncp) {
   y <- -ncp*x/sqrt(m+x^2)
   a <- (-y + sqrt(y^2 + 4*m)) / 2
   Hhmy <- 1/factorial(m) * a^m * exp(-0.5*(a+y)^2) * sqrt(2*pi*a^2/(m+a^2)) * (1 - 3*m/(4*(m+a^2)^2) + 5*m^2/(6*(m+a^2)^3))
factorial(m) / (2^((m-1)/2) * gamma(m/2) * sqrt(pi*m)) * exp(-0.5*m*ncp^2/(m+x^2)) * (m/(m+x^2))^((m+1)/2) * Hhmy
}

The function, as given, makes use of factorial() provided by the
gregmisc package (unless I overlooked it, I think it would be nice to
have a factorial function in the base package). The approximation part
is the calculation of Hhmy -- it is quite accurate though.

Actually, the two factorial(m) parts cancel each other out, so you can
rewrite the function as:

dtnc <- function(x, m, ncp) {
   y <- -ncp*x/sqrt(m+x^2)
   a <- (-y + sqrt(y^2 + 4*m)) / 2
   Hhmy <- a^m * exp(-0.5*(a+y)^2) * sqrt(2*pi*a^2/(m+a^2)) * (1 - 3*m/(4*(m+a^2)^2) + 5*m^2/(6*(m+a^2)^3))
1 / (2^((m-1)/2) * gamma(m/2) * sqrt(pi*m)) * exp(-0.5*m*ncp^2/(m+x^2)) * (m/(m+x^2))^((m+1)/2) * Hhmy
}

which doesn't require factorial(). I gave the first version, because the
exact value of Hhmy requires the evaluation on an integral, namely

integrand <- function(v, y, m) { v^m / factorial(m) * exp(-0.5*(v+y)^2) }
Hhmy      <- integrate(integrand, lower=0, upper=Inf, y=y, m=m)$value

and then, the two factorial(m)'s don't cancel. However, when using this
integral in the function above, it doesn't give the correct results. I
can't figure out why, but maybe this is still of interest to somebody.

--
Wolfgang Viechtbauer

"Are you not thinking what I am not thinking?"


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._