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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._