home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume34 / newmat06 / part01 / example.cxx next >
C/C++ Source or Header  |  1992-12-06  |  10KB  |  276 lines
  1. //$$ example.cxx                             Example of use of matrix package
  2.  
  3. #define WANT_STREAM                  // include.h will get stream fns
  4. #define WANT_MATH                    // include.h will get math fns
  5.  
  6. #include "include.h"                 // include standard files
  7.  
  8. #include "newmatap.h"                // need matrix applications
  9.  
  10. Real t3(Real);                       // round to 3 decimal places
  11.  
  12.  
  13. // demonstration of matrix package on linear regression problem
  14.  
  15.  
  16. void test1(Real* y, Real* x1, Real* x2, int nobs, int npred)
  17. {
  18.    cout << "\n\nTest 1 - traditional\n";
  19.  
  20.    // traditional sum of squares and products method of calculation
  21.    // with subtraction of means
  22.  
  23.    // make matrix of predictor values
  24.    Matrix X(nobs,npred);
  25.  
  26.    // load x1 and x2 into X
  27.    //    [use << rather than = with submatrices and/or loading arrays]
  28.    X.Column(1) << x1;  X.Column(2) << x2;
  29.  
  30.    // vector of Y values
  31.    ColumnVector Y(nobs); Y << y;
  32.  
  33.    // make vector of 1s
  34.    ColumnVector Ones(nobs); Ones = 1.0;
  35.  
  36.    // calculate means (averages) of x1 and x2 [ .t() takes transpose]
  37.    RowVector M = Ones.t() * X / nobs;
  38.  
  39.    // and subtract means from x1 and x1
  40.    Matrix XC(nobs,npred);
  41.    XC = X - Ones * M;
  42.  
  43.    // do the same to Y [need "Real" to convert 1x1 matrix to scalar]
  44.    ColumnVector YC(nobs);
  45.    Real m = (Ones.t() * Y).AsScalar() / nobs;  YC = Y - Ones * m;
  46.  
  47.    // form sum of squares and product matrix
  48.    //    [use << rather than = for copying Matrix into SymmetricMatrix]
  49.    SymmetricMatrix SSQ; SSQ << XC.t() * XC;
  50.  
  51.    // calculate estimate
  52.    //    [bracket last two terms to force this multiplication first]
  53.    //    [ .i() means inverse, but inverse is not explicity calculated]
  54.    ColumnVector A = SSQ.i() * (XC.t() * YC);
  55.  
  56.    // calculate estimate of constant term
  57.    Real a = m - (M * A).AsScalar();
  58.  
  59.    // Get variances of estimates from diagonal elements of invoice of SSQ
  60.    //    [ we are taking inverse of SSQ; would have been better to use
  61.    //        CroutMatrix method - see documentation ]
  62.    Matrix ISSQ = SSQ.i(); DiagonalMatrix D; D << ISSQ;
  63.    ColumnVector V = D.AsColumn();
  64.    Real v = 1.0/nobs + (M * ISSQ * M.t()).AsScalar();
  65.                         // for calc variance const
  66.  
  67.    // Calculate fitted values and residuals
  68.    int npred1 = npred+1;
  69.    ColumnVector Fitted = X * A + a;
  70.    ColumnVector Residual = Y - Fitted;
  71.    Real ResVar = Residual.SumSquare() / (nobs-npred1);
  72.  
  73.    // Get diagonals of Hat matrix (an expensive way of doing this)
  74.    Matrix X1(nobs,npred1); X1.Column(1)<<Ones; X1.Columns(2,npred1)<<X;
  75.    DiagonalMatrix Hat;  Hat << X1 * (X1.t() * X1).i() * X1.t();
  76.  
  77.    // print out answers
  78.    cout << "\nEstimates and their standard errors\n\n";
  79.    cout << a <<"\t"<< sqrt(v*ResVar) << "\n";
  80.    for (int i=1; i<=npred; i++)
  81.    cout << A(i) <<"\t"<< sqrt(V(i)*ResVar) << "\n";
  82.    cout << "\nObservations, fitted value, residual value, hat value\n";
  83.    for (i=1; i<=nobs; i++)
  84.       cout << X(i,1) <<"\t"<< X(i,2) <<"\t"<< Y(i) <<"\t"<<
  85.       t3(Fitted(i)) <<"\t"<< t3(Residual(i)) <<"\t"<< t3(Hat(i)) <<"\n";
  86.    cout << "\n\n";
  87. }
  88.  
  89. void test2(Real* y, Real* x1, Real* x2, int nobs, int npred)
  90. {
  91.    cout << "\n\nTest 2 - Cholesky\n";
  92.  
  93.    // traditional sum of squares and products method of calculation
  94.    // with subtraction of means - using Cholesky decomposition
  95.  
  96.    Matrix X(nobs,npred);
  97.    X.Column(1) << x1;  X.Column(2) << x2;
  98.    ColumnVector Y(nobs); Y << y;
  99.    ColumnVector Ones(nobs); Ones = 1.0;
  100.    RowVector M = Ones.t() * X / nobs;
  101.    Matrix XC(nobs,npred);
  102.    XC = X - Ones * M;
  103.    ColumnVector YC(nobs);
  104.    Real m = (Ones.t() * Y).AsScalar() / nobs;  YC = Y - Ones * m;
  105.    SymmetricMatrix SSQ; SSQ << XC.t() * XC;
  106.  
  107.    // Cholesky decomposition of SSQ
  108.    LowerTriangularMatrix L = Cholesky(SSQ);
  109.  
  110.    // calculate estimate
  111.    ColumnVector A = L.t().i() * (L.i() * (XC.t() * YC));
  112.  
  113.    // calculate estimate of constant term
  114.    Real a = m - (M * A).AsScalar();
  115.  
  116.    // Get variances of estimates from diagonal elements of invoice of SSQ
  117.    DiagonalMatrix D; D << L.t().i() * L.i();
  118.    ColumnVector V = D.AsColumn();
  119.    Real v = 1.0/nobs + (L.i() * M.t()).SumSquare();
  120.  
  121.    // Calculate fitted values and residuals
  122.    int npred1 = npred+1;
  123.    ColumnVector Fitted = X * A + a;
  124.    ColumnVector Residual = Y - Fitted;
  125.    Real ResVar = Residual.SumSquare() / (nobs-npred1);
  126.  
  127.    // Get diagonals of Hat matrix (an expensive way of doing this)
  128.    Matrix X1(nobs,npred1); X1.Column(1)<<Ones; X1.Columns(2,npred1)<<X;
  129.    DiagonalMatrix Hat;  Hat << X1 * (X1.t() * X1).i() * X1.t();
  130.  
  131.    // print out answers
  132.    cout << "\nEstimates and their standard errors\n\n";
  133.    cout << a <<"\t"<< sqrt(v*ResVar) << "\n";
  134.    for (int i=1; i<=npred; i++)
  135.       cout << A(i) <<"\t"<< sqrt(V(i)*ResVar) << "\n";
  136.    cout << "\nObservations, fitted value, residual value, hat value\n";
  137.    for (i=1; i<=nobs; i++)
  138.       cout << X(i,1) <<"\t"<< X(i,2) <<"\t"<< Y(i) <<"\t"<<
  139.       t3(Fitted(i)) <<"\t"<< t3(Residual(i)) <<"\t"<< t3(Hat(i)) <<"\n";
  140.    cout << "\n\n";
  141. }
  142.  
  143. void test3(Real* y, Real* x1, Real* x2, int nobs, int npred)
  144. {
  145.    cout << "\n\nTest 3 - Householder triangularisation\n";
  146.  
  147.    // Householder triangularisation method
  148.  
  149.    // load data - 1s into col 1 of matrix
  150.    int npred1 = npred+1;
  151.    Matrix X(nobs,npred1); ColumnVector Y(nobs);
  152.    X.Column(1) << 1.0;  X.Column(2) << x1;  X.Column(3) << x2;  Y << y;
  153.  
  154.    // do Householder triangularisation
  155.    // no need to deal with constant term separately
  156.    Matrix XT = X.t();             // Want data to be along rows
  157.    RowVector YT = Y.t();
  158.    LowerTriangularMatrix L; RowVector M;
  159.    HHDecompose(XT, L); HHDecompose(XT, YT, M); // YT now contains resids
  160.    ColumnVector A = L.t().i() * M.t();
  161.    ColumnVector Fitted = X * A;
  162.    Real ResVar = YT.SumSquare() / (nobs-npred1);
  163.  
  164.    // get variances of estimates
  165.    L = L.i(); DiagonalMatrix D; D << L.t() * L;
  166.  
  167.    // Get diagonals of Hat matrix
  168.    DiagonalMatrix Hat;  Hat << XT.t() * XT;
  169.  
  170.    // print out answers
  171.    cout << "\nEstimates and their standard errors\n\n";
  172.    for (int i=1; i<=npred1; i++)
  173.       cout << A(i) <<"\t"<< sqrt(D(i)*ResVar) << "\n";
  174.    cout << "\nObservations, fitted value, residual value, hat value\n";
  175.    for (i=1; i<=nobs; i++)
  176.       cout << X(i,2) <<"\t"<< X(i,3) <<"\t"<< Y(i) <<"\t"<<
  177.       t3(Fitted(i)) <<"\t"<< t3(YT(i)) <<"\t"<< t3(Hat(i)) <<"\n";
  178.    cout << "\n\n";
  179. }
  180.  
  181. void test4(Real* y, Real* x1, Real* x2, int nobs, int npred)
  182. {
  183.    cout << "\n\nTest 4 - singular value\n";
  184.  
  185.    // Singular value decomposition method
  186.  
  187.    // load data - 1s into col 1 of matrix
  188.    int npred1 = npred+1;
  189.    Matrix X(nobs,npred1); ColumnVector Y(nobs);
  190.    X.Column(1) << 1.0;  X.Column(2) << x1;  X.Column(3) << x2;  Y << y;
  191.  
  192.    // do SVD
  193.    Matrix U, V; DiagonalMatrix D;
  194.    SVD(X,D,U,V);                              // X = U * D * V.t()
  195.    ColumnVector Fitted = U.t() * Y;
  196.    ColumnVector A = V * ( D.i() * Fitted );
  197.    Fitted = U * Fitted;
  198.    ColumnVector Residual = Y - Fitted;
  199.    Real ResVar = Residual.SumSquare() / (nobs-npred1);
  200.  
  201.    // get variances of estimates
  202.    D << V * (D * D).i() * V.t();
  203.  
  204.    // Get diagonals of Hat matrix
  205.    DiagonalMatrix Hat;  Hat << U * U.t();
  206.  
  207.    // print out answers
  208.    cout << "\nEstimates and their standard errors\n\n";
  209.    for (int i=1; i<=npred1; i++)
  210.       cout << A(i) <<"\t"<< sqrt(D(i)*ResVar) << "\n";
  211.    cout << "\nObservations, fitted value, residual value, hat value\n";
  212.    for (i=1; i<=nobs; i++)
  213.       cout << X(i,2) <<"\t"<< X(i,3) <<"\t"<< Y(i) <<"\t"<<
  214.       t3(Fitted(i)) <<"\t"<< t3(Residual(i)) <<"\t"<< t3(Hat(i)) <<"\n";
  215.    cout << "\n\n";
  216. }
  217.  
  218. main()
  219. {
  220.    cout << "\nDemonstration of Matrix package\n\n";
  221.  
  222.    // Test for any memory not deallocated after running this program
  223.    Real* s1; { ColumnVector A(8000); s1 = A.Store(); }
  224.  
  225.    {
  226.       // the data
  227.       // you may need to read this data using cin if you are using a
  228.       // compiler that doesn't understand aggregates
  229. #ifndef ATandT
  230.       Real y[]  = { 8.3, 5.5, 8.0, 8.5, 5.7, 4.4, 6.3, 7.9, 9.1 };
  231.       Real x1[] = { 2.4, 1.8, 2.4, 3.0, 2.0, 1.2, 2.0, 2.7, 3.6 };
  232.       Real x2[] = { 1.7, 0.9, 1.6, 1.9, 0.5, 0.6, 1.1, 1.0, 0.5 };
  233. #else
  234.       Real y[9], x1[9], x2[9];
  235.       y[0]=8.3; y[1]=5.5; y[2]=8.0; y[3]=8.5; y[4]=5.7;
  236.       y[5]=4.4; y[6]=6.3; y[7]=7.9; y[8]=9.1;
  237.       x1[0]=2.4; x1[1]=1.8; x1[2]=2.4; x1[3]=3.0; x1[4]=2.0;
  238.       x1[5]=1.2; x1[6]=2.0; x1[7]=2.7; x1[8]=3.6;
  239.       x2[0]=1.7; x2[1]=0.9; x2[2]=1.6; x2[3]=1.9; x2[4]=0.5;
  240.       x2[5]=0.6; x2[6]=1.1; x2[7]=1.0; x2[8]=0.5;
  241. #endif
  242.       int nobs = 9;                           // number of observations
  243.       int npred = 2;                          // number of predictor values
  244.  
  245.       // we want to find the values of a,a1,a2 to give the best
  246.       // fit of y[i] with a0 + a1*x1[i] + a2*x2[i]
  247.       // Also print diagonal elements of hat matrix, X*(X.t()*X).i()*X.t()
  248.  
  249.       // this example demonstrates four methods of calculation
  250.  
  251.       Try
  252.       {
  253.          test1(y, x1, x2, nobs, npred);
  254.          test2(y, x1, x2, nobs, npred);
  255.          test3(y, x1, x2, nobs, npred);
  256.          test4(y, x1, x2, nobs, npred);
  257.       }
  258.       Catch(DataException) { cout << "\nInvalid data\n"; }
  259.       Catch(SpaceException) { cout << "\nMemory exhausted\n"; }
  260.       CatchAll { cout << "\nUnexpected program failure\n"; }
  261.    }
  262.  
  263. #ifdef DO_FREE_CHECK
  264.    FreeCheck::Status();
  265. #endif
  266.    Real* s2; { ColumnVector A(8000); s2 = A.Store(); }
  267.    cout << "\n\nChecking for lost memory: "
  268.       << (unsigned long)s1 << " " << (unsigned long)s2 << " ";
  269.    if (s1 != s2) cout << " - error\n"; else cout << " - ok\n";
  270.  
  271.    return 0;
  272.  
  273. }
  274.  
  275. Real t3(Real r) { return int(r*1000) / 1000.0; }
  276.