# [R] Accuracy problem in dchisq for non-central chi-squared

Peter Dalgaard BSA p.dalgaard at biostat.ku.dk
Thu Dec 14 11:56:34 CET 2000

```Uffe Høgsbro Thygesen <uht at dfu.min.dk> writes:

> Hi,
>
> I think I have identified a inaccuracy in dchisq when the non-centrality
> parameter is non-zero and large. Here's a little test:
>
> sys.dchisq.test <- function(N = 100000,mean = 0)
> {
>   z <- rnorm(N,mean = mean, sd = 1)
>   x <- z^2
>   xmin <- min(x)
>   xmax <- max(x)
>   br <- seq(xmin,xmax,length = 101)
>   dbr <- br[2]-br[1]
>   hist(x,br)
>   p <- dchisq(br,df = 1,ncp = mean^2) * dbr
>   lines(br,N*p)
> }
>
> par(mfrow = c(2,1))
> sys.dchisq.test(mean = 10)
> sys.dchisq.test(mean = 15)
>
> Here's my version information:
>
> platform Windows
> arch     x86
> os       Win32
> system   x86, Win32
> status
> major    1
> minor    1.1
> year     2000
> month    August
> day      15
> language R

Hej Uffe,

Same thing happens with 1.2pre, i.e. it is *not* fixed along with the
noncentral F problem. From the looks of the formula, I suspect that
the series expansion is getting terminated prematurely. Things work if
you use the definition directly; this works (slowly) up to at least
mean = 50:

function(N = 100000,mean = 0)
{
xnchisq<-function(x,df,lambda){N<-3*lambda;w<-dpois(0:N,lambda/2);
sapply(x,function(x)sum(w*dchisq(x,df=1+2*(0:N))))}
z <- rnorm(N,mean = mean, sd = 1)
x <- z^2
xmin <- min(x)
xmax <- max(x)
br <- seq(xmin,xmax,length = 101)
dbr <- br[2]-br[1]
hist(x,br)
p <- xnchisq(br,1,lambda = mean^2) * dbr
lines(br,N*p)
}

(BTW, overlaying histograms and desities is easier if you use freq=F
on the hist() call)
--
O__  ---- Peter Dalgaard             Blegdamsvej 3
c/ /'_ --- Dept. of Biostatistics     2200 Cph. N
(*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

```