MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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



  • Prev by Date: Re: integration for CDF[NoncentralChiSquareDistribution]
  • Next by Date: Re: bug in Integrate
  • Previous by thread: Re: integration for CDF[NoncentralChiSquareDistribution]
  • Next by thread: Re: Solveand Eliminatechoke on simple systems of equations