Home
Reading
Searching
Subscribe
Sponsors
Statistics
Posting
Contact
Spam
Lists
Links
About
Hosting
Filtering
Features Download
Marketing
Archives
FAQ
Blog
 
Gmane
From: <josef.pktd <at> gmail.com>
Subject: Re: Strange discontinuity in noncentral chisquare
Newsgroups: gmane.comp.python.scientific.user
Date: Thursday 28th May 2009 18:53:31 UTC (over 7 years ago)
On Thu, May 28, 2009 at 1:56 PM, Neal Becker  wrote:
> def pmiss2 (x, esnodB, N):
>    esno = 10**(0.1 * esnodB) * N
>    var = 1/esno
>    _lambda = 1/(0.5*var)
>
>    return ncx2.cdf (x, 2, _lambda)
>
> x = np.arange (0, 50, 0.1)
> p1 = [pmiss2 (e, 3.5, 24) for e in x]
>
> What's with this strange discontinuity?
> print p1:
> ...
> 3.475382846574262e-21,
>  4.2226227741362447e-21,
>  5.1248671653198949e-21,
>  6.2130949241675783e-21,
>  7.5242411687161146e-21,
>  9.1022970542215721e-21,
>  5.8787514615651664e-09,
>  6.2565279721932619e-09,
>  6.656924144742753e-09,
>  7.0811937641411923e-09,
>  7.5306544171300622e-09,
>  8.0066904400027596e-09,
>  8.5107559880675483e-09,
>  9.044378231221098e-09,
>  9.6091606801490766e-09,
> ...
>

This must be a bug in scipy.special.chndtr

class ncx2_gen(rv_continuous):

    def _cdf(self, x, df, nc):
        return special.chndtr(x,df,nc)

here is the isolated example

>>> ncx2.cdf (np.arange (20, 25, 0.2), 2, 1.07458615e+02)
array([  8.53614872e-23,   1.31445107e-22,   2.01359832e-22,
         3.06896021e-22,   4.65417198e-22,   7.02373115e-22,
         1.05489191e-21,   1.57689175e-21,   2.34632278e-21,
         3.47538283e-21,   5.12486714e-21,   7.52424113e-21,
         5.87875088e-09,   6.65692349e-09,   7.53065368e-09,
         8.51075516e-09,   9.60915975e-09,   1.08390218e-08,
         1.22148313e-08,   1.37525362e-08,   1.54696748e-08,
         1.73854823e-08,   1.95211903e-08,   2.18999761e-08,
         2.45472869e-08])

this uses numerical integration of the pdf, which is slow but doesn't
have the discontinuity

>>> ncx2.veccdf(np.arange (20, 25, 0.2), 2, 1.07458615e+02)
array([  1.21805117e-09,   1.39734318e-09,   1.60114784e-09,
         1.83255976e-09,   2.09503169e-09,   2.39241204e-09,
         2.72898607e-09,   3.10952087e-09,   3.53931444e-09,
         4.02424929e-09,   4.57085090e-09,   5.18635130e-09,
         5.87875834e-09,   6.65693102e-09,   7.53066129e-09,
         8.51076283e-09,   9.60916749e-09,   1.08390296e-08,
         1.22148392e-08,   1.37525442e-08,   1.54696828e-08,
         1.73855268e-08,   1.95212353e-08,   2.19000216e-08,
         2.45473328e-08])


Thanks for reporting, the test suite is not so fine tuned to catch these
cases.

The quick fix would be to replace the call to special with the
numerical integration, but this will make the cdf much slower ( for
the cases where it is correct). I did this for some other
distributions.
But for the use as the distribution of a test statistic the jump is
irrelevant, if you reject a hypotheses with 1e-9 or 1e-21 doesn't
really make a difference.



Josef
 
CD: 3ms