Name: David Williams E-mail Address: d.williams@qub.ac.uk IDL version: 5.4 Platform and OS: All Description of Behavior: The following is from an IDL newsgroup article by David on 16 November 2001.Example Code: None. Know Workarounds or Fixes: This workaround was suggested by Wayne Landsman on the same day David published the article above.This is a slightly more specialised problem than the usual general gripes, and it really has a big impact on our group's work, so I'm hoping someone else will have come across the same problem and found a remedy for it.
The code I use depends on the CHISQR_CVF and CHISQR_PDF functions to calculate confidence (or significance) levels for data that we analyse with wavelet transforms. These transforms aren't part of the IDL Wavelet Toolkit, because we've been running our own routines from before the time when wavelets were included in IDL. (If it's any help, the basis of the code is the WAVELET package written by Torrence & Compo, which you can download from http://paos.colorado.edu/research/wavelets .)
Our data often consist of thousands of points, and for any more than a 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.
Has anyone else (particularly those running similar wavelet analysis routines) run into these difficulties? And/or does anyone have a suggested solution or work-around, short of re-writing the CHISQR_PDF routine?
RSI Technical Support Response: They have been notified. No response as yet. Notes: Craig Markquart offers this as a note: Subject: Re: Problems with CHISQR_PDF when moving from IDL 5.2 to 5.4 From: Craig Markwardt Date: 16 Nov 2001 14:34:18 -0600 If you want to go ahead and take Wayne's suggestion, that's a good one. I've also just put up a set of routines for hypothesis testing and confidence interval estimation on my web page. This includes chi-square hypothesis testing. One major difference is that I took my code from Stephen Moshier's CEPHES special function library, so there's no dependence on the implementation of special functions by IDL. So far I haven't had convergence problems with this code. Good luck, Craig http://cow.physics.wisc.edu/~craigm/idl/ (under fitting) -------------------------------------------------------------------------- Craig B. Markwardt, Ph.D. EMAIL: craigmnet@cow.physics.wisc.edu --------------------------------------------------------------------------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.
whereasIDL> 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.