Re: integration for CDF[NoncentralChiSquareDistribution] failing

*To*: mathgroup at smc.vnet.net*Subject*: [mg74836] Re: integration for CDF[NoncentralChiSquareDistribution] failing*From*: "Ray Koopman" <koopman at sfu.ca>*Date*: Sat, 7 Apr 2007 04:05:17 -0400 (EDT)*References*: <ev4vp4$8ri$1@smc.vnet.net>

On Apr 6, 1:18 am, Albert Maydeu-Olivares <amay... at ub.edu> wrote: > Hi, > > The built-in function appears to work well for small number of > degrees of freedom. For large df, integration converges too slowly > and often fails. For instance > > 1 - CDF[NoncentralChiSquareDistribution[57, 5], 1] > > yields 1.04, when as a probability it should be < 1. Does anybody > have a function that works better than the built-in? Also, has > anybody implemented a normal approximation? I had planned to work > with 300+ df. > > Albert Here is some code for the cdf (F) and complementary cdf (G). It is based on the relation X[n,L] = X[n-1,0] + Z^2, where X[n,L] is a chi-square variable with df = n and noncentrality = L, and Z is normal with mean = Sqrt[L] and variance = 1. Then F[n,L,x] = Integrate[F[n-1,0,x-z^2]*f[z-Sqrt[L]], {z,-Sqrt[x],Sqrt[x]}], where f is the standard normal pdf. Although this is often slower than integrating the pdf, it's never very slow. L=0 and n=1 are special cases. I wrote the first version of this many years ago, and have tweaked it a few times since, but it's never been tested properly. Caveat user. F[n_,0 ,x_] := GammaRegularized[n/2,0,.5x]; F[n_,0.,x_] := GammaRegularized[n/2,0,.5x]; F[1 ,L_,x_] := .5 Erf[Sqrt[.5L]-Sqrt[.5x], Sqrt[.5L]+Sqrt[.5x]]; F[1.,L_,x_] := .5 Erf[Sqrt[.5L]-Sqrt[.5x], Sqrt[.5L]+Sqrt[.5x]]; F[n_,L_,x_] := With[{m = (n-1)/2, y = .5x, u = Sqrt[.5L], v = -Sqrt[.5L]}, NIntegrate[ GammaRegularized[m,0,(#*UnitStep@#&)[y-s^2]] * (Exp[-(s+u)^2]+Exp[-(s+v)^2]), {s,0,Sqrt[y]} ] / N@Sqrt@Pi ] G[n_,0 ,x_] := GammaRegularized[n/2,.5x]; G[n_,0.,x_] := GammaRegularized[n/2,.5x]; G[1 ,L_,x_] := .5 (Erfc[Sqrt[.5x]-Sqrt[.5L]] + Erfc[Sqrt[.5x]+Sqrt[.5L]]); G[1.,L_,x_] := .5 (Erfc[Sqrt[.5x]-Sqrt[.5L]] + Erfc[Sqrt[.5x]+Sqrt[.5L]]); G[n_,L_,x_] := With[{m = (n-1)/2, y = .5x, u = Sqrt[.5L], v = -Sqrt[.5L]}, NIntegrate[ GammaRegularized[m,(#*UnitStep@#&)[y-s^2]] * (Exp[-(s+u)^2]+Exp[-(s+v)^2]), {s,0,Sqrt[y]} ] / N@Sqrt@Pi + G[1,L,x] ]