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] ]