example.cpp
Go to the documentation of this file.
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                                      // newmatap.h will get include.h
00014 
00015 #include "newmatap.h"                // need matrix applications
00016 
00017 #include "newmatio.h"                // need matrix output routines
00018 
00019 #ifdef use_namespace
00020 using namespace NEWMAT;              // access NEWMAT namespace
00021 #endif
00022 
00023 
00024 // demonstration of matrix package on linear regression problem
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    // traditional sum of squares and products method of calculation
00032    // but not adjusting means; maybe subject to round-off error
00033 
00034    // make matrix of predictor values with 1s into col 1 of matrix
00035    int npred1 = npred+1;        // number of cols including col of ones.
00036    Matrix X(nobs,npred1);
00037    X.column(1) = 1.0;
00038 
00039    // load x1 and x2 into X
00040    //    [use << rather than = when loading arrays]
00041    X.column(2) << x1;  X.column(3) << x2;
00042 
00043    // vector of Y values
00044    ColumnVector Y(nobs); Y << y;
00045 
00046    // form sum of squares and product matrix
00047    //    [use << rather than = for copying Matrix into SymmetricMatrix]
00048    SymmetricMatrix SSQ; SSQ << X.t() * X;
00049 
00050    // calculate estimate
00051    //    [bracket last two terms to force this multiplication first]
00052    //    [ .i() means inverse, but inverse is not explicity calculated]
00053    ColumnVector A = SSQ.i() * (X.t() * Y);
00054 
00055    // Get variances of estimates from diagonal elements of inverse of SSQ
00056    // get inverse of SSQ - we need it for finding D
00057    DiagonalMatrix D; D << SSQ.i();
00058    ColumnVector V = D.as_column();
00059 
00060    // Calculate fitted values and residuals
00061    ColumnVector Fitted = X * A;
00062    ColumnVector Residual = Y - Fitted;
00063    Real ResVar = sum_square(Residual) / (nobs-npred1);
00064 
00065    // Get diagonals of Hat matrix (an expensive way of doing this)
00066    DiagonalMatrix Hat;  Hat << X * (X.t() * X).i() * X.t();
00067 
00068    // print out answers
00069    cout << "\nEstimates and their standard errors\n\n";
00070 
00071    // make vector of standard errors
00072    ColumnVector SE(npred1);
00073    for (int i=1; i<=npred1; i++) SE(i) = sqrt(V(i)*ResVar);
00074    // use concatenation function to form matrix and use matrix print
00075    // to get two columns
00076    cout << setw(11) << setprecision(5) << (A | SE) << endl;
00077 
00078    cout << "\nObservations, fitted value, residual value, hat value\n";
00079 
00080    // use concatenation again; select only columns 2 to 3 of X
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    // traditional sum of squares and products method of calculation
00091    // with subtraction of means - less subject to round-off error
00092    // than test1
00093 
00094    // make matrix of predictor values
00095    Matrix X(nobs,npred);
00096 
00097    // load x1 and x2 into X
00098    //    [use << rather than = when loading arrays]
00099    X.column(1) << x1;  X.column(2) << x2;
00100 
00101    // vector of Y values
00102    ColumnVector Y(nobs); Y << y;
00103 
00104    // make vector of 1s
00105    ColumnVector Ones(nobs); Ones = 1.0;
00106 
00107    // calculate means (averages) of x1 and x2 [ .t() takes transpose]
00108    RowVector M = X.sum_columns() / nobs;
00109 
00110    // and subtract means from x1 and x2
00111    Matrix XC(nobs,npred);
00112    XC = X - Ones * M;
00113 
00114    // do the same to Y [use sum to get sum of elements]
00115    ColumnVector YC(nobs);
00116    Real m = sum(Y) / nobs;  YC = Y - Ones * m;
00117 
00118    // form sum of squares and product matrix
00119    //    [use << rather than = for copying Matrix into SymmetricMatrix]
00120    SymmetricMatrix SSQ; SSQ << XC.t() * XC;
00121 
00122    // calculate estimate
00123    //    [bracket last two terms to force this multiplication first]
00124    //    [ .i() means inverse, but inverse is not explicity calculated]
00125    ColumnVector A = SSQ.i() * (XC.t() * YC);
00126 
00127    // calculate estimate of constant term
00128    //    [AsScalar converts 1x1 matrix to Real]
00129    Real a = m - (M * A).as_scalar();
00130 
00131    // Get variances of estimates from diagonal elements of inverse of SSQ
00132    //    [ we are taking inverse of SSQ - we need it for finding D ]
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                                             // for calc variance of const
00137 
00138    // Calculate fitted values and residuals
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    // Get diagonals of Hat matrix (an expensive way of doing this)
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    // print out answers
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    // make vector of standard errors
00154    ColumnVector SE(npred);
00155    for (int i=1; i<=npred; i++) SE(i) = sqrt(V(i)*ResVar);
00156    // use concatenation function to form matrix and use matrix print
00157    // to get two columns
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    // traditional sum of squares and products method of calculation
00170    // with subtraction of means - using Cholesky decomposition
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    // Cholesky decomposition of SSQ
00184    LowerTriangularMatrix L = Cholesky(SSQ);
00185 
00186    // calculate estimate
00187    ColumnVector A = L.t().i() * (L.i() * (XC.t() * YC));
00188 
00189    // calculate estimate of constant term
00190    Real a = m - (M * A).AsScalar();
00191 
00192    // Get variances of estimates from diagonal elements of invoice of SSQ
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    // Calculate fitted values and residuals
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    // Get diagonals of Hat matrix (an expensive way of doing this)
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    // print out answers
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    // QR triangularisation method
00226  
00227    // load data - 1s into col 1 of matrix
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    // do Householder triangularisation
00233    // no need to deal with constant term separately
00234    Matrix X1 = X;                 // Want copy of matrix
00235    ColumnVector Y1 = Y;
00236    UpperTriangularMatrix U; ColumnVector M;
00237    QRZ(X1, U); QRZ(X1, Y1, M);    // Y1 now contains resids
00238    ColumnVector A = U.i() * M;
00239    ColumnVector Fitted = X * A;
00240    Real ResVar = sum_square(Y1) / (nobs-npred1);
00241 
00242    // get variances of estimates
00243    U = U.i(); DiagonalMatrix D; D << U * U.t();
00244 
00245    // Get diagonals of Hat matrix
00246    DiagonalMatrix Hat;  Hat << X1 * X1.t();
00247 
00248    // print out answers
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    // Singular value decomposition method
00264  
00265    // load data - 1s into col 1 of matrix
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    // do SVD
00271    Matrix U, V; DiagonalMatrix D;
00272    SVD(X,D,U,V);                              // X = U * D * V.t()
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    // get variances of estimates
00280    D << V * (D * D).i() * V.t();
00281 
00282    // Get diagonals of Hat matrix
00283    DiagonalMatrix Hat;  Hat << U * U.t();
00284 
00285    // print out answers
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    // Test for any memory not deallocated after running this program
00302    Real* s1; { ColumnVector A(8000); s1 = A.data(); }
00303 
00304    {
00305       // the data
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;                           // number of observations
00313       int npred = 2;                          // number of predictor values
00314 
00315       // we want to find the values of a,a1,a2 to give the best
00316       // fit of y[i] with a0 + a1*x1[i] + a2*x2[i]
00317       // Also print diagonal elements of hat matrix, X*(X.t()*X).i()*X.t()
00318 
00319       // this example demonstrates five methods of calculation
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 


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:12