commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Phil Steitz <>
Subject Re: [math] Improving numerics in OLSMultipleLinearRegression
Date Tue, 10 Jun 2008 01:34:47 GMT
Mauro Talevi wrote:
> Hi Phil,
> thanks for reviewing the multiple linear regression implementations 
> and setting up the R/NIST data tests.   I finally got around to 
> installing R and can now run them too.
> Phil Steitz wrote:
>> While clear and elegant from a matrix algebra standpoint, the 
>> "nailve" implementation in OLSMultipleLinearRegression has bad 
>> numerical qualities.  It is well known that solving the normal 
>> equations directly does not give good numerics.  I just added some 
>> tests to actually verify parameter values, using the classic "Longly" 
>> dataset, for which NIST provides certified statistics.  This is a 
>> "hard" design matrix.  R was able to get to within 1E-8 of the 
>> certified parameter values.  OLSMultipleLinearRegression can only get 
>> 1E-1.
> The OLS implementation has been added as a simple by-product of the 
> GLS case - which is the main one I have needed for hypothesis testing 
> - as it came "for free" with unitary covariance.
> True - the emphasis was on clarity and formulaic simplicity.  And also 
> following the old Donald Knuth maxim "optimization is the root of all 
> evil".  But it seems like there is a need for refinement of the 
> implementation - the devil raised his head :-)
Yes, and I would distinguish performance optimization from numerical 
accuracy.  From my perspective, we can release a ".0" with room for 
performance improvement, but at least decent numerics are required.
>> We have talked in the past about providing an implementation based on 
>> QR decomposition.   Anyone up for  using the QR decomposition that we 
>> now have to do this?  I really think we need to do it (or something 
>> else to improve numerics) before releasing this class.  I will get to 
>> it eventually, but am a little pegged at the moment.  I will review 
>> and apply patches if someone is willing to do the implementation.  I 
>> can also explain here or offline how the R tests and NIST datasets 
>> work, as these are useful in validating code.
> I'd be happy to improve the impl.  I'm getting my head around R and 
> NIST, but perhaps a chat offline would not hurt!
I may be hard to catch synchronously, as my day-and-night job is a 
little demanding, but I would be happy to answer questions (with maybe a 
little latency ;)
>> Another thing that we should think about before releasing any of this 
>> stuff is the completeness of the API.  Many standard regression 
>> statistics are missing.  If we are going to stick with the Interface 
>> / Implementation setup, we need to get the right stuff into the 
>> interface.  It is also awkward to have to insert "1"'s in the design 
>> matrix to get an intercept term computed.  This is convenient for 
>> implementation, but awkward for users.  A more natural setup (IMHO) 
>> would be to expose a "noIntercept" or "hasIntercept" property for the 
>> model.
> No problem with adding other statistics - let's just decide on what is 
> the stardard regression API.
Here are some initial ideas on what should be included in the multiple 
regression API.  Other suggestions welcome!

1.  Coefficients should be accompanied by standard errors, t-statistics, 
two-sided t probablilities (can get these using t distribution from 
distributions package) and ideally confidence intervals.
2.  F, R-square, adjusted R-square, F prob (again can use distributions 
package to estimate)
3.  ANOVA table (Regression sum of squares, residual sum of squares)
4.  Residuals

R, SAS, SPSS and Excel all represent (or in the case of R, can 
construct) these basic statistics in some way in their output.  We 
should model them in classes representing properties of the computed 
> And finally, how do you see the no/hasIntercept model working?
As a configurable property - noIntercept means the model is estimated 
without an intercept.   The point I was making was more how the data is 
supplied via the API.  It is awkward to have to fill in a column of 1's 
to get the linear algebra to work to estimate a model with intercept 
(which should be the default).

I would recommend that we have setData or "newData" provide a n x m 
matrix, where n is the number of observations and m-1 is the number of 
independent variables.  Then either a) have the constructor take another 
argument specifying which column holds the dependent variable b) assume 
it is the first column c) support column labels and some form of model 
specification such as what R provides (a lot of work) d) split off the y 
vector, so setting data requires separate x and y vectors.  Probably a) 
is easiest for users, who will most often be starting with a rectangular 
array of data with the dependent variable in one of the columns.


To unsubscribe, e-mail:
For additional commands, e-mail:

View raw message