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

[Peter Dalgaard BSA]

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