R help - Dec 2000 - Accuracy problem in dchisq for non-central chi-squared

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         


Any comments on this one?

Best regards,

Uffe

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Uffe H. Thygesen
M.Sc.&Eng., Ph.D.
Danish Institute of Fisheries Research
http://www.dfu.min.dk


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

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

Uffe H?gsbro Thygesen <uht at dfu.min.dk> writes:
> Any comments on this one?
Ahem. I went to look at the code and it stared me right in the face:

    const static int maxiter = 100;

That needs to be keyed to the actual value of lambda! At lambda=225,
we're not even at the mode of the weighting Poisson distibution which
has mean lambda/2.

I think we can sneak a fix for that one in before the release of 1.2.0
tomorrow.  (pnchisq.c has 

    const static int itrmax = 10000;

which I should probably clone for now, but there must be a better
long-term solution for both of them. Anyone need a small exercise for
their students?)

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