[R] Non-central distributions

[Wolfgang Viechtbauer]

> all have a slot for the non-centrality parameter "ncp", of
> the functions for the t and F distributions:
>
>      dt(x, df, log = FALSE)
>      pt(q, df, ncp=0, lower.tail = TRUE, log.p = FALSE)
>      qt(p, df,        lower.tail = TRUE, log.p = FALSE)
>      rt(n, df)

Here is a function for calculating the density of a non-central t.

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

This is an approximation based on Resnikoff & Lieberman (1957). But it's
quite accurate. The exact pdf can either be expressed as an infinite sum
or requires the evaluation of an integral. I tried to implement the
latter using integrate(), but the results were very erratic. I checked
the accuracy of the approximation in various ways:

Comparing dtnc to dt in case ncp = 0 results in a maximum error of .0046
at the point x = 0 and df = 1. For df = 10, that error is already down
to .000021. I also compared pt with

integrate(dtnc, lower=-Inf, upper=x, df=df, ncp=ncp)

for a wide range of df's and ncp's and found that the results closely
matched to at least three decimal places even for small df. For 12 df,
differences only showed up in the 5th decimal. Obviously, use of
integrate() might contribute a little bit to the discrepancies.

Also, I checked whether dtnc integrates to 1 with

integrate(dtnc, lower=-Inf, upper=Inf, df=df, ncp=ncp)

for a wide range of df's and ncp's. Differences from 1 only showed up in
the second decimal places for df = 1, third decimal place for df = 2,
and the fourth decimal place for df = 4, and so on ...

So, it looks like the approximation provides pretty accurate results.

--
Wolfgang Viechtbauer

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