00001
00002
00003
00009
00010
00011 #define WANT_STREAM // include.h will get stream fns
00012 #define WANT_MATH // include.h will get math fns
00013
00014
00015 #include "newmatap.h"
00016
00017 #include "newmatio.h"
00018
00019 #ifdef use_namespace
00020 using namespace NEWMAT;
00021 #endif
00022
00023
00024
00025
00026
00027 void test1(Real* y, Real* x1, Real* x2, int nobs, int npred)
00028 {
00029 cout << "\n\nTest 1 - traditional, bad\n";
00030
00031
00032
00033
00034
00035 int npred1 = npred+1;
00036 Matrix X(nobs,npred1);
00037 X.column(1) = 1.0;
00038
00039
00040
00041 X.column(2) << x1; X.column(3) << x2;
00042
00043
00044 ColumnVector Y(nobs); Y << y;
00045
00046
00047
00048 SymmetricMatrix SSQ; SSQ << X.t() * X;
00049
00050
00051
00052
00053 ColumnVector A = SSQ.i() * (X.t() * Y);
00054
00055
00056
00057 DiagonalMatrix D; D << SSQ.i();
00058 ColumnVector V = D.as_column();
00059
00060
00061 ColumnVector Fitted = X * A;
00062 ColumnVector Residual = Y - Fitted;
00063 Real ResVar = sum_square(Residual) / (nobs-npred1);
00064
00065
00066 DiagonalMatrix Hat; Hat << X * (X.t() * X).i() * X.t();
00067
00068
00069 cout << "\nEstimates and their standard errors\n\n";
00070
00071
00072 ColumnVector SE(npred1);
00073 for (int i=1; i<=npred1; i++) SE(i) = sqrt(V(i)*ResVar);
00074
00075
00076 cout << setw(11) << setprecision(5) << (A | SE) << endl;
00077
00078 cout << "\nObservations, fitted value, residual value, hat value\n";
00079
00080
00081 cout << setw(9) << setprecision(3) <<
00082 (X.columns(2,3) | Y | Fitted | Residual | Hat.as_column());
00083 cout << "\n\n";
00084 }
00085
00086 void test2(Real* y, Real* x1, Real* x2, int nobs, int npred)
00087 {
00088 cout << "\n\nTest 2 - traditional, OK\n";
00089
00090
00091
00092
00093
00094
00095 Matrix X(nobs,npred);
00096
00097
00098
00099 X.column(1) << x1; X.column(2) << x2;
00100
00101
00102 ColumnVector Y(nobs); Y << y;
00103
00104
00105 ColumnVector Ones(nobs); Ones = 1.0;
00106
00107
00108 RowVector M = X.sum_columns() / nobs;
00109
00110
00111 Matrix XC(nobs,npred);
00112 XC = X - Ones * M;
00113
00114
00115 ColumnVector YC(nobs);
00116 Real m = sum(Y) / nobs; YC = Y - Ones * m;
00117
00118
00119
00120 SymmetricMatrix SSQ; SSQ << XC.t() * XC;
00121
00122
00123
00124
00125 ColumnVector A = SSQ.i() * (XC.t() * YC);
00126
00127
00128
00129 Real a = m - (M * A).as_scalar();
00130
00131
00132
00133 Matrix ISSQ = SSQ.i(); DiagonalMatrix D; D << ISSQ;
00134 ColumnVector V = D.AsColumn();
00135 Real v = 1.0/nobs + (M * ISSQ * M.t()).as_scalar();
00136
00137
00138
00139 int npred1 = npred+1;
00140 ColumnVector Fitted = X * A + a;
00141 ColumnVector Residual = Y - Fitted;
00142 Real ResVar = sum_square(Residual) / (nobs-npred1);
00143
00144
00145 Matrix X1(nobs,npred1); X1.column(1)<<Ones; X1.columns(2,npred1)<<X;
00146 DiagonalMatrix Hat; Hat << X1 * (X1.t() * X1).i() * X1.t();
00147
00148
00149 cout << "\nEstimates and their standard errors\n\n";
00150 cout.setf(ios::fixed, ios::floatfield);
00151 cout << setw(11) << setprecision(5) << a << " ";
00152 cout << setw(11) << setprecision(5) << sqrt(v*ResVar) << endl;
00153
00154 ColumnVector SE(npred);
00155 for (int i=1; i<=npred; i++) SE(i) = sqrt(V(i)*ResVar);
00156
00157
00158 cout << setw(11) << setprecision(5) << (A | SE) << endl;
00159 cout << "\nObservations, fitted value, residual value, hat value\n";
00160 cout << setw(9) << setprecision(3) <<
00161 (X | Y | Fitted | Residual | Hat.as_column());
00162 cout << "\n\n";
00163 }
00164
00165 void test3(Real* y, Real* x1, Real* x2, int nobs, int npred)
00166 {
00167 cout << "\n\nTest 3 - Cholesky\n";
00168
00169
00170
00171
00172 Matrix X(nobs,npred);
00173 X.column(1) << x1; X.column(2) << x2;
00174 ColumnVector Y(nobs); Y << y;
00175 RowVector M = X.sum_columns() / nobs;
00176 Matrix XC(nobs,npred);
00177 ColumnVector Ones(nobs); Ones = 1.0;
00178 XC = X - Ones * M;
00179 ColumnVector YC(nobs);
00180 Real m = sum(Y) / nobs; YC = Y - Ones * m;
00181 SymmetricMatrix SSQ; SSQ << XC.t() * XC;
00182
00183
00184 LowerTriangularMatrix L = Cholesky(SSQ);
00185
00186
00187 ColumnVector A = L.t().i() * (L.i() * (XC.t() * YC));
00188
00189
00190 Real a = m - (M * A).AsScalar();
00191
00192
00193 DiagonalMatrix D; D << L.t().i() * L.i();
00194 ColumnVector V = D.as_column();
00195 Real v = 1.0/nobs + sum_square(L.i() * M.t());
00196
00197
00198 int npred1 = npred+1;
00199 ColumnVector Fitted = X * A + a;
00200 ColumnVector Residual = Y - Fitted;
00201 Real ResVar = sum_square(Residual) / (nobs-npred1);
00202
00203
00204 Matrix X1(nobs,npred1); X1.column(1)<<Ones; X1.columns(2,npred1)<<X;
00205 DiagonalMatrix Hat; Hat << X1 * (X1.t() * X1).i() * X1.t();
00206
00207
00208 cout << "\nEstimates and their standard errors\n\n";
00209 cout.setf(ios::fixed, ios::floatfield);
00210 cout << setw(11) << setprecision(5) << a << " ";
00211 cout << setw(11) << setprecision(5) << sqrt(v*ResVar) << endl;
00212 ColumnVector SE(npred);
00213 for (int i=1; i<=npred; i++) SE(i) = sqrt(V(i)*ResVar);
00214 cout << setw(11) << setprecision(5) << (A | SE) << endl;
00215 cout << "\nObservations, fitted value, residual value, hat value\n";
00216 cout << setw(9) << setprecision(3) <<
00217 (X | Y | Fitted | Residual | Hat.as_column());
00218 cout << "\n\n";
00219 }
00220
00221 void test4(Real* y, Real* x1, Real* x2, int nobs, int npred)
00222 {
00223 cout << "\n\nTest 4 - QR triangularisation\n";
00224
00225
00226
00227
00228 int npred1 = npred+1;
00229 Matrix X(nobs,npred1); ColumnVector Y(nobs);
00230 X.column(1) = 1.0; X.column(2) << x1; X.column(3) << x2; Y << y;
00231
00232
00233
00234 Matrix X1 = X;
00235 ColumnVector Y1 = Y;
00236 UpperTriangularMatrix U; ColumnVector M;
00237 QRZ(X1, U); QRZ(X1, Y1, M);
00238 ColumnVector A = U.i() * M;
00239 ColumnVector Fitted = X * A;
00240 Real ResVar = sum_square(Y1) / (nobs-npred1);
00241
00242
00243 U = U.i(); DiagonalMatrix D; D << U * U.t();
00244
00245
00246 DiagonalMatrix Hat; Hat << X1 * X1.t();
00247
00248
00249 cout << "\nEstimates and their standard errors\n\n";
00250 ColumnVector SE(npred1);
00251 for (int i=1; i<=npred1; i++) SE(i) = sqrt(D(i)*ResVar);
00252 cout << setw(11) << setprecision(5) << (A | SE) << endl;
00253 cout << "\nObservations, fitted value, residual value, hat value\n";
00254 cout << setw(9) << setprecision(3) <<
00255 (X.columns(2,3) | Y | Fitted | Y1 | Hat.as_column());
00256 cout << "\n\n";
00257 }
00258
00259 void test5(Real* y, Real* x1, Real* x2, int nobs, int npred)
00260 {
00261 cout << "\n\nTest 5 - singular value\n";
00262
00263
00264
00265
00266 int npred1 = npred+1;
00267 Matrix X(nobs,npred1); ColumnVector Y(nobs);
00268 X.column(1) = 1.0; X.column(2) << x1; X.column(3) << x2; Y << y;
00269
00270
00271 Matrix U, V; DiagonalMatrix D;
00272 SVD(X,D,U,V);
00273 ColumnVector Fitted = U.t() * Y;
00274 ColumnVector A = V * ( D.i() * Fitted );
00275 Fitted = U * Fitted;
00276 ColumnVector Residual = Y - Fitted;
00277 Real ResVar = sum_square(Residual) / (nobs-npred1);
00278
00279
00280 D << V * (D * D).i() * V.t();
00281
00282
00283 DiagonalMatrix Hat; Hat << U * U.t();
00284
00285
00286 cout << "\nEstimates and their standard errors\n\n";
00287 ColumnVector SE(npred1);
00288 for (int i=1; i<=npred1; i++) SE(i) = sqrt(D(i)*ResVar);
00289 cout << setw(11) << setprecision(5) << (A | SE) << endl;
00290 cout << "\nObservations, fitted value, residual value, hat value\n";
00291 cout << setw(9) << setprecision(3) <<
00292 (X.columns(2,3) | Y | Fitted | Residual | Hat.as_column());
00293 cout << "\n\n";
00294 }
00295
00296 int main()
00297 {
00298 cout << "\nDemonstration of Matrix package\n";
00299 cout << "\nPrint a real number (may help lost memory test): " << 3.14159265 << "\n";
00300
00301
00302 Real* s1; { ColumnVector A(8000); s1 = A.data(); }
00303
00304 {
00305
00306
00307 Real y[] = { 8.3, 5.5, 8.0, 8.5, 5.7, 4.4, 6.3, 7.9, 9.1 };
00308 Real x1[] = { 2.4, 1.8, 2.4, 3.0, 2.0, 1.2, 2.0, 2.7, 3.6 };
00309 Real x2[] = { 1.7, 0.9, 1.6, 1.9, 0.5, 0.6, 1.1, 1.0, 0.5 };
00310
00311
00312 int nobs = 9;
00313 int npred = 2;
00314
00315
00316
00317
00318
00319
00320
00321 Try
00322 {
00323 test1(y, x1, x2, nobs, npred);
00324 test2(y, x1, x2, nobs, npred);
00325 test3(y, x1, x2, nobs, npred);
00326 test4(y, x1, x2, nobs, npred);
00327 test5(y, x1, x2, nobs, npred);
00328 }
00329 CatchAll { cout << BaseException::what(); }
00330 }
00331
00332 #ifdef DO_FREE_CHECK
00333 FreeCheck::Status();
00334 #endif
00335 Real* s2; { ColumnVector A(8000); s2 = A.data(); }
00336 cout << "\n\nThe following test does not work with all compilers - see documentation\n";
00337 cout << "Checking for lost memory: "
00338 << (unsigned long)s1 << " " << (unsigned long)s2 << " ";
00339 if (s1 != s2) cout << " - see section 2.8\n"; else cout << " - ok\n";
00340
00341 return 0;
00342
00343 }
00344
00345