Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: integration for CDF[NoncentralChiSquareDistribution]

• To: mathgroup at smc.vnet.net
• Subject: [mg74831] Re: [mg74806] integration for CDF[NoncentralChiSquareDistribution]
• From: Darren Glosemeyer <darreng at wolfram.com>
• Date: Sat, 7 Apr 2007 04:02:44 -0400 (EDT)
• References: <200704060817.EAA09037@smc.vnet.net>

Albert Maydeu-Olivares 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
>   

Hi Albert,
This is an issue that has been addressed in the version of Mathematica 
currently under development. CDF relies on numeric integration of the 
pdf for NoncentralChiSquareDistribution (note that in the version under 
development at least one of the numeric values will need to be inexact 
to make CDF use numeric methods). In effect, the source of the problem 
was numeric instability introduced by some auto-expansion of the pdf.

The solution is to only evaluate the pdf for numeric values, thus 
avoiding the symbolic expansion. In version 5.2, the results can be 
obtained as follows.


In[1]:= << Statistics`

In[2]:= pdffun[nn_?NumericQ, ll_?NumericQ,t_?NumericQ] :=
        (E^(-ll/2 - t/2)* (t/ll)^((-2 + nn)/4)*BesselI[(-2 + nn)/2, 
Sqrt[ll*t]])/2;


This gives the CDF at 1.

In[3]:= NIntegrate[pdffun[57, 5, x], {x, 0, 1}]

                 -41
Out[3]= 8.5132 10


Note that with a machine precision result this value is too small to be 
noticed in the subtraction from 1.

In[4]:= 1-%

Out[4]= 1.


Switching to a higher precision computation, the difference from 1 can 
be seen.

In[5]:= NIntegrate[pdffun[57, 5, x], {x, 0, 1},WorkingPrecision->20]

                      -41
Out[5]= 8.513197603 10

In[6]:= 1-%

Out[6]= 0.99999999999999999999999999999999999999991486802397


This will also allow for computation for large numbers of degrees of 
freedom. However, if you are interested in a normal approximation, 
Johnson, Kotz and Balakrishnan give several normal approximations in 
Continuous Univariate Distributions, Volume 2, and you could use one of 
those which is appropriate for the degrees of freedom and noncentrality 
values in your application.

Darren Glosemeyer
Wolfram Research

• References:
  • integration for CDF[NoncentralChiSquareDistribution] failing
    • From: Albert Maydeu-Olivares <amaydeu@ub.edu>