[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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list