Java Code Examples for org.apache.commons.math3.linear.MatrixUtils#createRealIdentityMatrix()
The following examples show how to use
org.apache.commons.math3.linear.MatrixUtils#createRealIdentityMatrix() .
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
| Source File: GLSMultipleLinearRegressionTest.java From astor with GNU General Public License v2.0 | 6 votes |
|
/**
* Verifies that GLS with identity covariance matrix gives the same results
* as OLS.
*/
@Test
public void testGLSOLSConsistency() {
RealMatrix identityCov = MatrixUtils.createRealIdentityMatrix(16);
GLSMultipleLinearRegression glsModel = new GLSMultipleLinearRegression();
OLSMultipleLinearRegression olsModel = new OLSMultipleLinearRegression();
glsModel.newSampleData(longley, 16, 6);
olsModel.newSampleData(longley, 16, 6);
glsModel.newCovarianceData(identityCov.getData());
double[] olsBeta = olsModel.calculateBeta().toArray();
double[] glsBeta = glsModel.calculateBeta().toArray();
// TODO: Should have assertRelativelyEquals(double[], double[], eps) in TestUtils
// Should also add RealVector and RealMatrix versions
for (int i = 0; i < olsBeta.length; i++) {
TestUtils.assertRelativelyEquals(olsBeta[i], glsBeta[i], 10E-7);
}
}
Example 2
| Source File: GLSMultipleLinearRegressionTest.java From astor with GNU General Public License v2.0 | 6 votes |
|
/**
* Verifies that GLS with identity covariance matrix gives the same results
* as OLS.
*/
@Test
public void testGLSOLSConsistency() {
RealMatrix identityCov = MatrixUtils.createRealIdentityMatrix(16);
GLSMultipleLinearRegression glsModel = new GLSMultipleLinearRegression();
OLSMultipleLinearRegression olsModel = new OLSMultipleLinearRegression();
glsModel.newSampleData(longley, 16, 6);
olsModel.newSampleData(longley, 16, 6);
glsModel.newCovarianceData(identityCov.getData());
double[] olsBeta = olsModel.calculateBeta().toArray();
double[] glsBeta = glsModel.calculateBeta().toArray();
// TODO: Should have assertRelativelyEquals(double[], double[], eps) in TestUtils
// Should also add RealVector and RealMatrix versions
for (int i = 0; i < olsBeta.length; i++) {
TestUtils.assertRelativelyEquals(olsBeta[i], glsBeta[i], 10E-7);
}
}
Example 3
| Source File: GLSMultipleLinearRegressionTest.java From astor with GNU General Public License v2.0 | 6 votes |
|
/**
* Verifies that GLS with identity covariance matrix gives the same results
* as OLS.
*/
@Test
public void testGLSOLSConsistency() {
RealMatrix identityCov = MatrixUtils.createRealIdentityMatrix(16);
GLSMultipleLinearRegression glsModel = new GLSMultipleLinearRegression();
OLSMultipleLinearRegression olsModel = new OLSMultipleLinearRegression();
glsModel.newSampleData(longley, 16, 6);
olsModel.newSampleData(longley, 16, 6);
glsModel.newCovarianceData(identityCov.getData());
double[] olsBeta = olsModel.calculateBeta().toArray();
double[] glsBeta = glsModel.calculateBeta().toArray();
// TODO: Should have assertRelativelyEquals(double[], double[], eps) in TestUtils
// Should also add RealVector and RealMatrix versions
for (int i = 0; i < olsBeta.length; i++) {
TestUtils.assertRelativelyEquals(olsBeta[i], glsBeta[i], 10E-7);
}
}
Example 4
| Source File: GLSMultipleLinearRegressionTest.java From astor with GNU General Public License v2.0 | 6 votes |
|
/**
* Verifies that GLS with identity covariance matrix gives the same results
* as OLS.
*/
@Test
public void testGLSOLSConsistency() {
RealMatrix identityCov = MatrixUtils.createRealIdentityMatrix(16);
GLSMultipleLinearRegression glsModel = new GLSMultipleLinearRegression();
OLSMultipleLinearRegression olsModel = new OLSMultipleLinearRegression();
glsModel.newSampleData(longley, 16, 6);
olsModel.newSampleData(longley, 16, 6);
glsModel.newCovarianceData(identityCov.getData());
double[] olsBeta = olsModel.calculateBeta().toArray();
double[] glsBeta = glsModel.calculateBeta().toArray();
// TODO: Should have assertRelativelyEquals(double[], double[], eps) in TestUtils
// Should also add RealVector and RealMatrix versions
for (int i = 0; i < olsBeta.length; i++) {
TestUtils.assertRelativelyEquals(olsBeta[i], glsBeta[i], 10E-7);
}
}
Example 5
| Source File: IntegerCopyNumberTransitionProbabilityCacheCollectionUnitTest.java From gatk-protected with BSD 3-Clause "New" or "Revised" License | 5 votes |
|
private static RealMatrix getDirectMatrixPower(final RealMatrix mat, final int power) {
if (power < 0) {
throw new IllegalArgumentException("Can not calculate negative matrix powers");
} else if (power == 0) {
return MatrixUtils.createRealIdentityMatrix(mat.getColumnDimension());
} else {
return mat.power(power);
}
}
Example 6
| Source File: WeightedIntDiGraphTest.java From pacaya with Apache License 2.0 | 5 votes |
|
@Test
public void testSumWalks() {
RealMatrix m = new Array2DRowRealMatrix(new double[][] {
{ 0.3, 0.0, 0.3 },
{ 0.1, 0.6, 0.7 },
{ 0.0, 0.4, 0.2 },
});
RealMatrix sum = MatrixUtils.createRealIdentityMatrix(3);
for (int i = 0; i < 2000; i++) {
sum = sum.add(m.power(i + 1));
}
RealVector ones = new ArrayRealVector(3, 1.0);
// this is how I'm computing the reference (different from the methods implementation)
double expectedSumWalks = sum.preMultiply(ones).dotProduct(ones);
RealVector s = new ArrayRealVector(new double[] {1, 2, 3});
RealVector t = new ArrayRealVector(new double[] {5, 7, 4});
assertEquals(127.49999999999903, expectedSumWalks, tol);
assertEquals(127.49999999999903, WeightedIntDiGraph.sumWalks(m, null, null), tol);
assertEquals(127.49999999999903, WeightedIntDiGraph.sumWalks(m, null, ones), tol);
assertEquals(127.49999999999903, WeightedIntDiGraph.sumWalks(m, ones, null), tol);
assertTrue(TestUtils.checkThrows(() -> {
WeightedIntDiGraph.sumWalks(
new Array2DRowRealMatrix(
new double[][] {
{ 0.3, 0.0, 0.3, 0 },
{ 0.1, 0.6, 0.7, 0 },
{ 0.0, 0.4, 0.2, 0 }}),
null,
null);
}, InputMismatchException.class));
assertEquals(699.9999999999948, WeightedIntDiGraph.sumWalks(m, null, t), tol);
assertEquals(274.99999999999795, WeightedIntDiGraph.sumWalks(m, s, null), tol);
assertEquals(1509.9999999999886, WeightedIntDiGraph.sumWalks(m, s, t), tol);
}
Example 7
| Source File: KalmanFilter.java From astor with GNU General Public License v2.0 | 5 votes |
|
/**
* Correct the current state estimate with an actual measurement.
*
* @param z
* the measurement vector
* @throws NullArgumentException
* if the measurement vector is {@code null}
* @throws DimensionMismatchException
* if the dimension of the measurement vector does not fit
* @throws SingularMatrixException
* if the covariance matrix could not be inverted
*/
public void correct(final RealVector z)
throws NullArgumentException, DimensionMismatchException, SingularMatrixException {
// sanity checks
MathUtils.checkNotNull(z);
if (z.getDimension() != measurementMatrix.getRowDimension()) {
throw new DimensionMismatchException(z.getDimension(),
measurementMatrix.getRowDimension());
}
// S = H * P(k) - * H' + R
RealMatrix s = measurementMatrix.multiply(errorCovariance)
.multiply(measurementMatrixT)
.add(measurementModel.getMeasurementNoise());
// invert S
// as the error covariance matrix is a symmetric positive
// semi-definite matrix, we can use the cholesky decomposition
DecompositionSolver solver = new CholeskyDecomposition(s).getSolver();
RealMatrix invertedS = solver.getInverse();
// Inn = z(k) - H * xHat(k)-
RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation));
// calculate gain matrix
// K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1
// K(k) = P(k)- * H' * S^-1
RealMatrix kalmanGain = errorCovariance.multiply(measurementMatrixT).multiply(invertedS);
// update estimate with measurement z(k)
// xHat(k) = xHat(k)- + K * Inn
stateEstimation = stateEstimation.add(kalmanGain.operate(innovation));
// update covariance of prediction error
// P(k) = (I - K * H) * P(k)-
RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension());
errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance);
}
Example 8
| Source File: OLSMultipleLinearRegressionTest.java From astor with GNU General Public License v2.0 | 4 votes |
|
/**
* Test hat matrix computation
*
*/
@Test
public void testHat() {
/*
* This example is from "The Hat Matrix in Regression and ANOVA",
* David C. Hoaglin and Roy E. Welsch,
* The American Statistician, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
*
*/
double[] design = new double[] {
11.14, .499, 11.1,
12.74, .558, 8.9,
13.13, .604, 8.8,
11.51, .441, 8.9,
12.38, .550, 8.8,
12.60, .528, 9.9,
11.13, .418, 10.7,
11.7, .480, 10.5,
11.02, .406, 10.5,
11.41, .467, 10.7
};
int nobs = 10;
int nvars = 2;
// Estimate the model
OLSMultipleLinearRegression model = new OLSMultipleLinearRegression();
model.newSampleData(design, nobs, nvars);
RealMatrix hat = model.calculateHat();
// Reference data is upper half of symmetric hat matrix
double[] referenceData = new double[] {
.418, -.002, .079, -.274, -.046, .181, .128, .222, .050, .242,
.242, .292, .136, .243, .128, -.041, .033, -.035, .004,
.417, -.019, .273, .187, -.126, .044, -.153, .004,
.604, .197, -.038, .168, -.022, .275, -.028,
.252, .111, -.030, .019, -.010, -.010,
.148, .042, .117, .012, .111,
.262, .145, .277, .174,
.154, .120, .168,
.315, .148,
.187
};
// Check against reference data and verify symmetry
int k = 0;
for (int i = 0; i < 10; i++) {
for (int j = i; j < 10; j++) {
Assert.assertEquals(referenceData[k], hat.getEntry(i, j), 10e-3);
Assert.assertEquals(hat.getEntry(i, j), hat.getEntry(j, i), 10e-12);
k++;
}
}
/*
* Verify that residuals computed using the hat matrix are close to
* what we get from direct computation, i.e. r = (I - H) y
*/
double[] residuals = model.estimateResiduals();
RealMatrix I = MatrixUtils.createRealIdentityMatrix(10);
double[] hatResiduals = I.subtract(hat).operate(model.getY()).toArray();
TestUtils.assertEquals(residuals, hatResiduals, 10e-12);
}
Example 9
| Source File: BiDiagonalTransformerTest.java From astor with GNU General Public License v2.0 | 4 votes |
|
private void checkOrthogonal(RealMatrix m) {
RealMatrix mTm = m.transpose().multiply(m);
RealMatrix id = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
Assert.assertEquals(0, mTm.subtract(id).getNorm(), 1.0e-14);
}
Example 10
| Source File: OLSMultipleLinearRegressionTest.java From astor with GNU General Public License v2.0 | 4 votes |
|
/**
* Test hat matrix computation
*
*/
@Test
public void testHat() {
/*
* This example is from "The Hat Matrix in Regression and ANOVA",
* David C. Hoaglin and Roy E. Welsch,
* The American Statistician, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
*
*/
double[] design = new double[] {
11.14, .499, 11.1,
12.74, .558, 8.9,
13.13, .604, 8.8,
11.51, .441, 8.9,
12.38, .550, 8.8,
12.60, .528, 9.9,
11.13, .418, 10.7,
11.7, .480, 10.5,
11.02, .406, 10.5,
11.41, .467, 10.7
};
int nobs = 10;
int nvars = 2;
// Estimate the model
OLSMultipleLinearRegression model = new OLSMultipleLinearRegression();
model.newSampleData(design, nobs, nvars);
RealMatrix hat = model.calculateHat();
// Reference data is upper half of symmetric hat matrix
double[] referenceData = new double[] {
.418, -.002, .079, -.274, -.046, .181, .128, .222, .050, .242,
.242, .292, .136, .243, .128, -.041, .033, -.035, .004,
.417, -.019, .273, .187, -.126, .044, -.153, .004,
.604, .197, -.038, .168, -.022, .275, -.028,
.252, .111, -.030, .019, -.010, -.010,
.148, .042, .117, .012, .111,
.262, .145, .277, .174,
.154, .120, .168,
.315, .148,
.187
};
// Check against reference data and verify symmetry
int k = 0;
for (int i = 0; i < 10; i++) {
for (int j = i; j < 10; j++) {
Assert.assertEquals(referenceData[k], hat.getEntry(i, j), 10e-3);
Assert.assertEquals(hat.getEntry(i, j), hat.getEntry(j, i), 10e-12);
k++;
}
}
/*
* Verify that residuals computed using the hat matrix are close to
* what we get from direct computation, i.e. r = (I - H) y
*/
double[] residuals = model.estimateResiduals();
RealMatrix I = MatrixUtils.createRealIdentityMatrix(10);
double[] hatResiduals = I.subtract(hat).operate(model.getY()).toArray();
TestUtils.assertEquals(residuals, hatResiduals, 10e-12);
}
Example 11
| Source File: BiDiagonalTransformerTest.java From astor with GNU General Public License v2.0 | 4 votes |
|
private void checkOrthogonal(RealMatrix m) {
RealMatrix mTm = m.transpose().multiply(m);
RealMatrix id = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
Assert.assertEquals(0, mTm.subtract(id).getNorm(), 1.0e-14);
}
Example 12
| Source File: OLSMultipleLinearRegressionTest.java From astor with GNU General Public License v2.0 | 4 votes |
|
/**
* Test hat matrix computation
*
*/
@Test
public void testHat() {
/*
* This example is from "The Hat Matrix in Regression and ANOVA",
* David C. Hoaglin and Roy E. Welsch,
* The American Statistician, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
*
*/
double[] design = new double[] {
11.14, .499, 11.1,
12.74, .558, 8.9,
13.13, .604, 8.8,
11.51, .441, 8.9,
12.38, .550, 8.8,
12.60, .528, 9.9,
11.13, .418, 10.7,
11.7, .480, 10.5,
11.02, .406, 10.5,
11.41, .467, 10.7
};
int nobs = 10;
int nvars = 2;
// Estimate the model
OLSMultipleLinearRegression model = new OLSMultipleLinearRegression();
model.newSampleData(design, nobs, nvars);
RealMatrix hat = model.calculateHat();
// Reference data is upper half of symmetric hat matrix
double[] referenceData = new double[] {
.418, -.002, .079, -.274, -.046, .181, .128, .222, .050, .242,
.242, .292, .136, .243, .128, -.041, .033, -.035, .004,
.417, -.019, .273, .187, -.126, .044, -.153, .004,
.604, .197, -.038, .168, -.022, .275, -.028,
.252, .111, -.030, .019, -.010, -.010,
.148, .042, .117, .012, .111,
.262, .145, .277, .174,
.154, .120, .168,
.315, .148,
.187
};
// Check against reference data and verify symmetry
int k = 0;
for (int i = 0; i < 10; i++) {
for (int j = i; j < 10; j++) {
Assert.assertEquals(referenceData[k], hat.getEntry(i, j), 10e-3);
Assert.assertEquals(hat.getEntry(i, j), hat.getEntry(j, i), 10e-12);
k++;
}
}
/*
* Verify that residuals computed using the hat matrix are close to
* what we get from direct computation, i.e. r = (I - H) y
*/
double[] residuals = model.estimateResiduals();
RealMatrix I = MatrixUtils.createRealIdentityMatrix(10);
double[] hatResiduals = I.subtract(hat).operate(model.getY()).toArray();
TestUtils.assertEquals(residuals, hatResiduals, 10e-12);
}
Example 13
| Source File: OLSMultipleLinearRegressionTest.java From astor with GNU General Public License v2.0 | 4 votes |
|
/**
* Test hat matrix computation
*
* @throws Exception
*/
@Test
public void testHat() throws Exception {
/*
* This example is from "The Hat Matrix in Regression and ANOVA",
* David C. Hoaglin and Roy E. Welsch,
* The American Statistician, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
*
*/
double[] design = new double[] {
11.14, .499, 11.1,
12.74, .558, 8.9,
13.13, .604, 8.8,
11.51, .441, 8.9,
12.38, .550, 8.8,
12.60, .528, 9.9,
11.13, .418, 10.7,
11.7, .480, 10.5,
11.02, .406, 10.5,
11.41, .467, 10.7
};
int nobs = 10;
int nvars = 2;
// Estimate the model
OLSMultipleLinearRegression model = new OLSMultipleLinearRegression();
model.newSampleData(design, nobs, nvars);
RealMatrix hat = model.calculateHat();
// Reference data is upper half of symmetric hat matrix
double[] referenceData = new double[] {
.418, -.002, .079, -.274, -.046, .181, .128, .222, .050, .242,
.242, .292, .136, .243, .128, -.041, .033, -.035, .004,
.417, -.019, .273, .187, -.126, .044, -.153, .004,
.604, .197, -.038, .168, -.022, .275, -.028,
.252, .111, -.030, .019, -.010, -.010,
.148, .042, .117, .012, .111,
.262, .145, .277, .174,
.154, .120, .168,
.315, .148,
.187
};
// Check against reference data and verify symmetry
int k = 0;
for (int i = 0; i < 10; i++) {
for (int j = i; j < 10; j++) {
Assert.assertEquals(referenceData[k], hat.getEntry(i, j), 10e-3);
Assert.assertEquals(hat.getEntry(i, j), hat.getEntry(j, i), 10e-12);
k++;
}
}
/*
* Verify that residuals computed using the hat matrix are close to
* what we get from direct computation, i.e. r = (I - H) y
*/
double[] residuals = model.estimateResiduals();
RealMatrix I = MatrixUtils.createRealIdentityMatrix(10);
double[] hatResiduals = I.subtract(hat).operate(model.getY()).toArray();
TestUtils.assertEquals(residuals, hatResiduals, 10e-12);
}
Example 14
| Source File: BiDiagonalTransformerTest.java From astor with GNU General Public License v2.0 | 4 votes |
|
private void checkOrthogonal(RealMatrix m) {
RealMatrix mTm = m.transpose().multiply(m);
RealMatrix id = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
Assert.assertEquals(0, mTm.subtract(id).getNorm(), 1.0e-14);
}
Example 15
| Source File: MatrixStack.java From NOVA-Core with GNU Lesser General Public License v3.0 | 4 votes |
|
/**
* Replaces current transformation matrix by an identity matrix.
*/
public void loadIdentity() {
current = MatrixUtils.createRealIdentityMatrix(4);
}
Example 16
| Source File: KalmanFilter.java From astor with GNU General Public License v2.0 | 4 votes |
|
/**
* Correct the current state estimate with an actual measurement.
*
* @param z
* the measurement vector
* @throws NullArgumentException
* if the measurement vector is {@code null}
* @throws DimensionMismatchException
* if the dimension of the measurement vector does not fit
* @throws SingularMatrixException
* if the covariance matrix could not be inverted
*/
public void correct(final RealVector z)
throws NullArgumentException, DimensionMismatchException, SingularMatrixException {
// sanity checks
MathUtils.checkNotNull(z);
if (z.getDimension() != measurementMatrix.getRowDimension()) {
throw new DimensionMismatchException(z.getDimension(),
measurementMatrix.getRowDimension());
}
// S = H * P(k) * H' + R
RealMatrix s = measurementMatrix.multiply(errorCovariance)
.multiply(measurementMatrixT)
.add(measurementModel.getMeasurementNoise());
// Inn = z(k) - H * xHat(k)-
RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation));
// calculate gain matrix
// K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1
// K(k) = P(k)- * H' * S^-1
// instead of calculating the inverse of S we can rearrange the formula,
// and then solve the linear equation A x X = B with A = S', X = K' and B = (H * P)'
// K(k) * S = P(k)- * H'
// S' * K(k)' = H * P(k)-'
RealMatrix kalmanGain = new CholeskyDecomposition(s).getSolver()
.solve(measurementMatrix.multiply(errorCovariance.transpose()))
.transpose();
// update estimate with measurement z(k)
// xHat(k) = xHat(k)- + K * Inn
stateEstimation = stateEstimation.add(kalmanGain.operate(innovation));
// update covariance of prediction error
// P(k) = (I - K * H) * P(k)-
RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension());
errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance);
}
Example 17
| Source File: BiDiagonalTransformerTest.java From astor with GNU General Public License v2.0 | 4 votes |
|
private void checkOrthogonal(RealMatrix m) {
RealMatrix mTm = m.transpose().multiply(m);
RealMatrix id = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
Assert.assertEquals(0, mTm.subtract(id).getNorm(), 1.0e-14);
}
Example 18
| Source File: OLSMultipleLinearRegressionTest.java From astor with GNU General Public License v2.0 | 4 votes |
|
/**
* Test hat matrix computation
*
*/
@Test
public void testHat() {
/*
* This example is from "The Hat Matrix in Regression and ANOVA",
* David C. Hoaglin and Roy E. Welsch,
* The American Statistician, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
*
*/
double[] design = new double[] {
11.14, .499, 11.1,
12.74, .558, 8.9,
13.13, .604, 8.8,
11.51, .441, 8.9,
12.38, .550, 8.8,
12.60, .528, 9.9,
11.13, .418, 10.7,
11.7, .480, 10.5,
11.02, .406, 10.5,
11.41, .467, 10.7
};
int nobs = 10;
int nvars = 2;
// Estimate the model
OLSMultipleLinearRegression model = new OLSMultipleLinearRegression();
model.newSampleData(design, nobs, nvars);
RealMatrix hat = model.calculateHat();
// Reference data is upper half of symmetric hat matrix
double[] referenceData = new double[] {
.418, -.002, .079, -.274, -.046, .181, .128, .222, .050, .242,
.242, .292, .136, .243, .128, -.041, .033, -.035, .004,
.417, -.019, .273, .187, -.126, .044, -.153, .004,
.604, .197, -.038, .168, -.022, .275, -.028,
.252, .111, -.030, .019, -.010, -.010,
.148, .042, .117, .012, .111,
.262, .145, .277, .174,
.154, .120, .168,
.315, .148,
.187
};
// Check against reference data and verify symmetry
int k = 0;
for (int i = 0; i < 10; i++) {
for (int j = i; j < 10; j++) {
Assert.assertEquals(referenceData[k], hat.getEntry(i, j), 10e-3);
Assert.assertEquals(hat.getEntry(i, j), hat.getEntry(j, i), 10e-12);
k++;
}
}
/*
* Verify that residuals computed using the hat matrix are close to
* what we get from direct computation, i.e. r = (I - H) y
*/
double[] residuals = model.estimateResiduals();
RealMatrix I = MatrixUtils.createRealIdentityMatrix(10);
double[] hatResiduals = I.subtract(hat).operate(model.getY()).toArray();
TestUtils.assertEquals(residuals, hatResiduals, 10e-12);
}
Example 19
| Source File: OLSMultipleLinearRegressionTest.java From astor with GNU General Public License v2.0 | 4 votes |
|
/**
* Test hat matrix computation
*
*/
@Test
public void testHat() {
/*
* This example is from "The Hat Matrix in Regression and ANOVA",
* David C. Hoaglin and Roy E. Welsch,
* The American Statistician, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
*
*/
double[] design = new double[] {
11.14, .499, 11.1,
12.74, .558, 8.9,
13.13, .604, 8.8,
11.51, .441, 8.9,
12.38, .550, 8.8,
12.60, .528, 9.9,
11.13, .418, 10.7,
11.7, .480, 10.5,
11.02, .406, 10.5,
11.41, .467, 10.7
};
int nobs = 10;
int nvars = 2;
// Estimate the model
OLSMultipleLinearRegression model = new OLSMultipleLinearRegression();
model.newSampleData(design, nobs, nvars);
RealMatrix hat = model.calculateHat();
// Reference data is upper half of symmetric hat matrix
double[] referenceData = new double[] {
.418, -.002, .079, -.274, -.046, .181, .128, .222, .050, .242,
.242, .292, .136, .243, .128, -.041, .033, -.035, .004,
.417, -.019, .273, .187, -.126, .044, -.153, .004,
.604, .197, -.038, .168, -.022, .275, -.028,
.252, .111, -.030, .019, -.010, -.010,
.148, .042, .117, .012, .111,
.262, .145, .277, .174,
.154, .120, .168,
.315, .148,
.187
};
// Check against reference data and verify symmetry
int k = 0;
for (int i = 0; i < 10; i++) {
for (int j = i; j < 10; j++) {
Assert.assertEquals(referenceData[k], hat.getEntry(i, j), 10e-3);
Assert.assertEquals(hat.getEntry(i, j), hat.getEntry(j, i), 10e-12);
k++;
}
}
/*
* Verify that residuals computed using the hat matrix are close to
* what we get from direct computation, i.e. r = (I - H) y
*/
double[] residuals = model.estimateResiduals();
RealMatrix I = MatrixUtils.createRealIdentityMatrix(10);
double[] hatResiduals = I.subtract(hat).operate(model.getY()).toArray();
TestUtils.assertEquals(residuals, hatResiduals, 10e-12);
}
