Phil,
Phil Steitz wrote:
I think R uses QR as described above. Comments or suggestions for other
default implementations are most welcome. We should aim to provide a
default implementation that is reasonably fast and provides good
numerics across a broad range of design matrices.
Got around to testing QR decomposition on OLS.
The short answer is that is does not seems to make much difference.
Rather it looks like the dataset and the number of observations that are
far more significant for good numerics, as is also shown by your recent
addition of Swiss Fertility dataset (with nobs 3 times as large), for
which results match (either with or without QR) up to 10^-12 tolerance.
Here's the resulting numerics for comparison:
Longley dataset (nobs=16):
PL=[-3482258.6569276676, 15.06187677821299, -0.03581918037047249,
-2.0202298136474104, -1.0332268695603801, -0.051104103746114404,
1829.1514737363977]
QR=[-3482258.7119702557, 15.061873615257795, -0.03581918168712586,
-2.020229840231328, -1.0332268778552742, -0.05110409751647271,
1829.1515061042903]
LG=[-3482258.63459582, 15.0618722713733, -0.035819179292591,
-2.02022980381683, -1.03322686717359, -0.0511041056535807, 1829.15146461355]
Swiss Fertility dataset (nobs=47):
PL=[91.05542390271336, -0.22064551045713723, -0.26058239824327045,
-0.9616123845602972, 0.12441843147162471]
QR=[91.05542390271366, -0.22064551045714642, -0.26058239824326457,
-0.9616123845602974, 0.12441843147162669]
SF=[91.05542390271397, -0.22064551045715, -0.26058239824328,
-0.9616123845603, 0.12441843147162]
(Legend: PL = plain OLS, QR = QR-decomposed OLS, LG = Longley R results,
SF = Swiss Fertility R results).
Interestingly, it's only on the intercepts (ie the first regression
parameter) that we get the very poor numerics. While not a numerical
argument, one could say that the statistically more significant
parameter is the slope.
Anyway, attached is patch with QR-based implementation and modified test
to print out comparison results.
Cheers
Index:
src/test/org/apache/commons/math/stat/regression/OLSMultipleLinearRegressionTest.java
===================================================================
---
src/test/org/apache/commons/math/stat/regression/OLSMultipleLinearRegressionTest.java
(revision 669849)
+++
src/test/org/apache/commons/math/stat/regression/OLSMultipleLinearRegressionTest.java
(working copy)
@@ -16,9 +16,11 @@
*/
package org.apache.commons.math.stat.regression;
+import java.util.Arrays;
+
+import org.apache.commons.math.TestUtils;
import org.junit.Before;
import org.junit.Test;
-import org.apache.commons.math.TestUtils;
public class OLSMultipleLinearRegressionTest extends
AbstractMultipleLinearRegressionTest {
@@ -109,7 +111,18 @@
// Estimate the model
MultipleLinearRegression model = new OLSMultipleLinearRegression();
model.addData(y, x, null);
+
+
System.out.println("OL="+Arrays.toString(model.estimateRegressionParameters()));
+
+ MultipleLinearRegression qr = new QROLSMultipleLinearRegression();
+ qr.addData(y, x, null);
+
+
System.out.println("QR="+Arrays.toString(qr.estimateRegressionParameters()));
+ System.out.println("LG="+Arrays.toString(new
double[]{-3482258.63459582, 15.0618722713733,
+ -0.358191792925910E-01,-2.02022980381683,
+ -1.03322686717359,-0.511041056535807E-01,
+ 1829.15146461355}));
// Check expected beta values from NIST
double[] betaHat = model.estimateRegressionParameters();
TestUtils.assertEquals(betaHat,
@@ -201,6 +214,19 @@
MultipleLinearRegression model = new OLSMultipleLinearRegression();
model.addData(y, x, null);
+
System.out.println("OL="+Arrays.toString(model.estimateRegressionParameters()));
+
+ MultipleLinearRegression qr = new QROLSMultipleLinearRegression();
+ qr.addData(y, x, null);
+
+
System.out.println("QR="+Arrays.toString(qr.estimateRegressionParameters()));
+
+ System.out.println("LG="+Arrays.toString(new
double[]{91.05542390271397,
+ -0.22064551045715,
+ -0.26058239824328,
+ -0.96161238456030,
+ 0.12441843147162}));
+
// Check expected beta values from R
double[] betaHat = model.estimateRegressionParameters();
TestUtils.assertEquals(betaHat,
Index:
src/java/org/apache/commons/math/stat/regression/QROLSMultipleLinearRegression.java
===================================================================
---
src/java/org/apache/commons/math/stat/regression/QROLSMultipleLinearRegression.java
(revision 0)
+++
src/java/org/apache/commons/math/stat/regression/QROLSMultipleLinearRegression.java
(revision 0)
@@ -0,0 +1,107 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.stat.regression;
+
+import org.apache.commons.math.linear.QRDecomposition;
+import org.apache.commons.math.linear.QRDecompositionImpl;
+import org.apache.commons.math.linear.RealMatrix;
+
+
+/**
+ * The OLS implementation of the multiple linear regression.
+ *
+ * OLS assumes the covariance matrix of the error to be diagonal and with
equal variance.
+ * <pre>
+ * u ~ N(0, sigma^2*I)
+ * </pre>
+ *
+ * By QR decomposition, X = QR, Q'Q = I
+ *
+ * <pre>
+ * z = Q' y = R b
+ * </pre>
+ *
+ * Estimated by OLS,
+ * <pre>
+ * b=(X'X)^-1X'y = (R'R)^-1R'z = (R'R)^-1R'Q'y
+ * </pre>
+ * whose variance is
+ * <pre>
+ * Var(b)=MSE*(X'X)^-1, MSE=u'u/(n-k)
+ * </pre>
+ * @version $Revision: 658645 $ $Date: 2008-05-21 13:13:27 +0100 (Wed, 21 May
2008) $
+ * @since 2.0
+ */
+public class QROLSMultipleLinearRegression extends
AbstractMultipleLinearRegression {
+
+ private RealMatrix Q = null;
+ private RealMatrix R = null;
+ private RealMatrix Z = null;
+
+ /**
+ * [EMAIL PROTECTED]
+ */
+ public void addData(double[] y, double[][] x, double[][] covariance) {
+ validateSampleData(x, y);
+ addYSampleData(y);
+ addXSampleData(x);
+ QRDecomposition qr = new QRDecompositionImpl(X);
+ Q = qr.getQ();
+ R = qr.getR();
+ Z = Q.transpose().multiply(Y);
+ }
+
+ /**
+ * Calculates beta by OLS, using QR decomposition of X
+ * <pre>
+ * b=(X'X)^-1X'y = (R'R)^-1R'z
+ * </pre>
+ * @return beta
+ */
+ protected RealMatrix calculateBeta() {
+ RealMatrix RTR = R.transpose().multiply(R);
+ RealMatrix beta = RTR.inverse().multiply(R.transpose()).multiply(Z);
+ return beta;
+ }
+
+ /**
+ * Calculates the variance on the beta by OLS, using QR decomposition of X
+ * <pre>
+ * Var(b)=(X'X)^-1=(R'R)^-1
+ * </pre>
+ * @return The beta variance
+ */
+ protected RealMatrix calculateBetaVariance() {
+ RealMatrix RTR = R.transpose().multiply(R);
+ return RTR.inverse();
+ }
+
+
+ /**
+ * Calculates the variance on the Y by OLS.
+ * <pre>
+ * Var(y)=Tr(u'u)/(n-k)
+ * </pre>
+ * @return The Y variance
+ */
+ protected double calculateYVariance() {
+ RealMatrix u = calculateResiduals();
+ RealMatrix sse = u.transpose().multiply(u);
+ return sse.getTrace()/(X.getRowDimension()-X.getColumnDimension());
+ }
+
+}
---------------------------------------------------------------------
To unsubscribe, e-mail: [EMAIL PROTECTED]
For additional commands, e-mail: [EMAIL PROTECTED]
