Volume 9, Issue 4
Tricks of the Trade
In and Out
Download This Issue
Staff and Contributors
On the Numerical Accuracy of Mathematica 5.0 for Doing Linear and Nonlinear Regression
There is a small but growing literature on the numerical accuracy of various computer programming packages. A major contributor to the assessment of the accuracy of econometric software, as well as statistical add-ons to third-level programming languages widely used by econometricians, such as GAUSS, MATLAB, Maple, and Mathematica, has been Bruce D. McCullough ( and the references cited therein).
McCullough and Wilson [2, 27] assessed the reliability of the statistical procedures in Excel 97 for estimation of both linear and nonlinear regression, random number generation, and the calculation of cumulative distributions and concluded that "Excel's performance in all three areas is found to be inadequate. Persons desiring to conduct analyses of statistical data are advised not to use Excel." McCullough made a complete assessment of all of Mathematica 4.0's then extant statistical functions . In 2000, I made a similar assessment for ordinary least squares (OLS) linear regression analysis using Excel 2000 and the Numerical Algorithms Group's (NAG's) Statistical Add-Ins for Excel software (1). I found that the "NAG-enhanced" version of Excel did not do significantly better than the non-enhanced version. In the same paper, I also compared the performance of Mathematica 4.1.
The basis for McCullough's assessments as well as mine were comparisons with the benchmark results for the Statistical Reference Datasets (StRD) that are available from the NIST (2).
My comparison of Excel 2000 with the benchmark linear regression analyses produced results that did not differ significantly from those McCullough obtained for Excel 97. Notwithstanding the excellent reputation of NAG's software, the performance of the Excel add-on is not superior to Excel. For cases in which Excel performs well, the NAG add-on also does well, but when Excel performs unsatisfactorily, the NAG add-on does poorly as well. However, often the results differ from the certified results in a different direction from those obtained directly in Excel. In the case of linear regression, on the other hand, Excel on its own does not do much worse than many other packages analyzed by McCullough and others. My comparisons at the time suggested that McCullough's judgment, at least with respect to linear regression, may have been unduly critical; however, I do agree that the statistical functions in Excel should not be used with or without NAG enhancement.
Vinod [4, 211] undertook a careful analysis of the statistical functions in GAUSS for Windows Version 3.2.27, also using comparisons with the StRD. He found "... that the algorithms used by GAUSS sometimes fail rather badly to provide accurate results. We blame both the algorithms used and the language itself." I have not yet undertaken a similar reassessment of the more recent and allegedly significantly improved versions of GAUSS.
McCullough [5, 152] applied similar methodology to SAS, SPSS, and S-Plus, which I call fourth-level programming languages or statistical packages. With respect to univariate statistics, analysis of variance, and linear regression, McCullough reported that "All [SAS] results are quite accurate," in terms of six digits of accuracy or better, although SAS failed to deliver any solution in one of the worst cases. Even though SAS failed to achieve convergence in a number of cases, the results for nonlinear regression were similar. The performance of all three packages with respect to random number generation and the calculation of probabilities for statistical distributions was notably poorer than for OLS linear regression.
McCullough  analyzed the accuracy of Mathematica 4.0's statistics add-on package (3). Mathematica is not restricted as to precision; however, precision may be specified. The default is to use machine arithmetic, which is typically 16 digits of precision in internal calculations (a typical value of $MachinePrecision is or approximately 16). Of course this does not guarantee 16 digits of precision in the final result. Note that accuracy and precision are not the same as defined in Mathematica. A higher level of precision may be attained by setting $MinPrecision above 16 or by setting WorkingPrecision->n, where n is a number greater than MachinePrecision. Only a small number of Mathematica functions support the WorkingPrecision option. Moreover, when software other than Mathematica is used, increasing the level of precision in this way may not be possible. This is because within the code the computation has been carried out in floating-point notation. No matter how many degrees of arbitrary precision are requested, the result is still contaminated with the errors from floating-point arithmetic. Arbitrary precision only works well when the computation has been carried out uncontaminated with floating-point error. In Mathematica it is possible to carry out the computation with numbers expressed in rational form, which yields exact results, but this still will not avoid the problem of inverting an ill-conditioned matrix.
At default precision, Mathematica is quite accurate although McCullough obtained no solution in one StRD case of OLS linear regression and one in an StRD case of nonlinear regression. When $MinPrecision is raised to 20 digits, McCullough reports that astonishing, almost perfect, accuracy, as compared with the NIST benchmark results, is obtained in all cases considered. Actually, Mathematica's default level of precision gives exact matches to the benchmark results in most cases, but the results Mathematica displays are rounded to six places and thus may not look quite the same. To display the results to more places, use NumberForm. The increase in computational time using a higher level of precision than default in direct, and therefore symbolically pure, computation of the OLS estimates appears to be negligible no matter how high the level of precision is specified. Although it is certainly true that increased accuracy may not matter for the conclusions to be drawn from a statistical analysis, it may also happen that rounding errors in the course of an otherwise straightforward calculation may result in wildly different results or that no solution is found. If the cost of obtaining accurate results even in problematic cases is not great, it is certainly prudent to obtain what increased accuracy may be possible. If a preprogrammed package is used, as it is in this article, it may not be possible to control for rounding error because of the original coding.
In this article, I discuss the details of my own comparisons of Mathematica 5.0's standard statistical add-on package. McCullough employs summary measures, which must of necessity leave out a lot of relevant detail in the case of many parameters, but I present a full range of results. A separate section deals with the more computationally difficult benchmark cases for nonlinear regression.
In the next section, I describe the NIST StRD benchmarks, particularly the 11 designed for benchmarking OLS linear regression (i.e., linear in the parameters to be estimated, so polynomial regression is not excluded) and the 27 nonlinear benchmark datasets. I also discuss the way in which the NIST certified results have been obtained for the linear regression StRD set. In Section 3, I discuss the general problems of nonlinear regression analysis and the difficulty obtaining appropriate significance tests short of computing the parameter estimates by maximum likelihood (ML). In Section 4, I summarize the results of using Version 5.0 of Mathematica's Statistics`LinearRegression`standard add-on package and compare these with the StRD benchmark results. Essentially the same results may be obtained directly by writing out the standard textbook formulae in Mathematica, although these formulae are not the best way of carrying out the requisite numerical calculations. Details of the numerical results are presented in Appendices A and B, in which tabular comparisons are made of the NIST StRD benchmark results with those obtained using Mathematica's Statistics`LinearRegression` standard add-on package and those using the NonlinearRegress option in the Statistics`NonlinearFit` standard add-on package. The appendices are part of the electronic version of this article, but are not included in the printed version.
About Mathematica | Download Mathematica Player
© 2005 Wolfram Media, Inc. All rights reserved.