[SciPy-user] Strange discontinuity in noncentral chisquare

Robert Kern robert.kern@gmail....
Thu May 28 14:03:25 CDT 2009


On Thu, May 28, 2009 at 13:53,  <josef.pktd@gmail.com> wrote:
> On Thu, May 28, 2009 at 1:56 PM, Neal Becker <ndbecker2@gmail.com> 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

I notice the following snippets of code, which appear guilty:

C     .. Statement Functions ..
      LOGICAL qsmall
C     ..
C     .. Statement Function definitions ..
      qsmall(xx) = sum .LT. 1.0D-20 .OR. xx .LT. eps*sum

# That is a feature of Fortran I knew nothing about.

   60 IF (qsmall(term)) GO TO 80

   80 cum = sum
      RETURN

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco


More information about the SciPy-user mailing list