Re: Problems with CHISQR_PDF when moving from IDL 5.2 to 5.4 [message #28037 is a reply to message #28035] |
Fri, 16 November 2001 09:34  |
Wayne Landsman
Messages: 117 Registered: January 1997
|
Senior Member |
|
|
David Williams wrote:
> few hundred degrees of freedom (and we have thousands of degrees) IDL
> 5.4 runs into difficulties. The problem we're encountering is that
> CHISQR_PDF can't return an answer, and the only major difference I can
> immediately see -- between the version of CHISQR_PDF bundled with 5.2
> and that bundled with 5.4 -- is that the newer version calls IGAMMA
> rather than IGAMMA_PDF. In fact, it's IGAMMA that fails to converge to a
> value and echoes a message saying this. It's pretty frustrating because
> 5.2 handled our routines just fine.
You've isolated the problem, though since IGAMMA takes scalar parameters,
you can narrow in on the problem even further by simply noting that in V5.4
IDL> print,igamma(1300,1252)
returns % IGAMMA: Failed to converge within given parameters.
whereas IDL> print,igamma_pdf(1300,1252) returns the correct value of
0.0903564
IGAMMA uses the same iterative algorithm as the old IGAMMA_PDF except that
the default number of iterations is 100 instead of 1000. So I would
either change the line in the V5.4 CHISQR_PDF that reads
gres = igamma(df/2.0, (x > 0)/2.0)
to
gres = igamma(df/2.0, (x > 0)/2.0, ITMAX = 1000)
or change the default value of ITMAX in IGAMMA() to 1000.
(You can use the V5.2 version of chisqr_pdf.pro, though in general, except
for the ITMAX parameter, IGAMMA() looks more robust than IGAMMA_PDF())
The number of required iterations increases with the number of degrees of
freedom ( which determines the "a" parameter send to the GAMMA function)
and so in principle one could have ITMAX scale with the number of degrees
of freedom. But this is a fast scalar calculation, (and the iterations
stop as soon as one achieves the desired tolerance) so I see no reason to
not to give ITMAX a large value.
Wayne Landsman landsman@mpb.gsfc.nasa.gov
|
|
|