home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume21 / newmat02 / part03 / example.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1991-08-01  |  6.1 KB  |  175 lines
  1. //$$ example.cxx                             Example of use of matrix package
  2.  
  3. typedef double real;                 // all floating point double
  4.  
  5. #include "include.hxx"               // include standard files
  6.  
  7. #include "newmatap.hxx"              // need matrix applications
  8.  
  9. real t3(real);                       // round to 3 decimal places
  10.  
  11. // demonstration of matrix package on linear regression problem
  12.  
  13. main()
  14. {
  15.    // the data
  16.    // you may need to read this data using cin if you are using a
  17.    // compiler that doesn't understand aggregates
  18.    real y[]  = { 8.3, 5.5, 8.0, 8.5, 5.7, 4.4, 6.3, 7.9, 9.1 };
  19.    real x1[] = { 2.4, 1.8, 2.4, 3.0, 2.0, 1.2, 2.0, 2.7, 3.6 };
  20.    real x2[] = { 1.7, 0.9, 1.6, 1.9, 0.5, 0.6, 1.1, 1.0, 0.5 };
  21.    int nobs = 9;                           // number of observations
  22.    int npred = 2;                          // number of predictor values
  23.    int npred1 = npred+1;
  24.  
  25.    // we want to find the values of a,a1,a2 to give the best
  26.    // fit of y[i] with a0 + a1*x1[i] + a2*x2[i]
  27.  
  28.    // this example demonstrates three methods of calculation
  29.  
  30.    {
  31.       // traditional sum of squares and products method of calculation
  32.       // with subtraction of means
  33.  
  34.       // make matrix of predictor values
  35.       Matrix X(nobs,npred);
  36.  
  37.       // load x1 and x2 into X
  38.       //    [use << rather than = with submatrices and/or loading arrays]
  39.       X.Col(1) << x1;  X.Col(2) << x2;
  40.  
  41.       // vector of Y values
  42.       ColumnVector Y(nobs); Y << y;
  43.  
  44.       // make vector of 1s
  45.       ColumnVector Ones(nobs); Ones = 1.0;
  46.  
  47.       // calculate means (averages) of x1 and x2 [ .t() takes transpose]
  48.       RowVector M = Ones.t() * X / nobs;
  49.  
  50.       // and subtract means from x1 and x1
  51.       Matrix XC(nobs,npred);
  52.       XC = X - Ones * M;
  53.  
  54.       // do the same to Y [need (1,1) to convert 1x1 matrix to scalar]
  55.       ColumnVector YC(nobs);
  56.       real m = Matrix(Ones.t() * Y)(1,1) / nobs;  YC = Y - Ones * m;
  57.  
  58.       // form sum of squares and product matrix
  59.       //    [use << rather than = for copying Matrix into SymmetricMatrix]
  60.       SymmetricMatrix SSQ; SSQ << XC.t() * XC;
  61.  
  62.       // calculate estimate
  63.       //    [bracket last two terms to force this multiplication first]
  64.       //    [ .i() means inverse, but inverse is not explicity calculated]
  65.       ColumnVector A = SSQ.i() * (XC.t() * YC);
  66.  
  67.       // calculate estimate of constant term
  68.       real a = m - Matrix(M * A)(1,1);
  69.  
  70.       // Get variances of estimates from diagonal elements of invoice of SSQ
  71.       //    [ we are taking inverse of SSQ; would have been better to use
  72.       //        CroutMatrix method - see documentation ]
  73.       Matrix ISSQ = SSQ.i(); DiagonalMatrix D; D << ISSQ;
  74.       ColumnVector V = D.CopyToColumn();
  75.       real v = 1.0/nobs + Matrix(M * ISSQ * M.t())(1,1);
  76.                         // for calc variance const
  77.  
  78.       // Calculate fitted values and residuals
  79.       ColumnVector Fitted = X * A + a;
  80.       ColumnVector Residual = Y - Fitted;
  81.       real ResVar = Residual.SumSquare() / (nobs-npred1);
  82.  
  83.       // print out answers
  84.       cout << "\nEstimates and their standard errors\n\n";
  85.       cout << a <<"\t"<< sqrt(v*ResVar) << "\n";
  86.       for (int i=1; i<=npred; i++)
  87.      cout << A(i) <<"\t"<< sqrt(V(i)*ResVar) << "\n";
  88.       cout << "\nObservations, fitted value, residual value\n";
  89.       for (i=1; i<=nobs; i++)
  90.      cout << X(i,1) <<"\t"<< X(i,2) <<"\t"<< Y(i) <<"\t"<<
  91.      t3(Fitted(i)) <<"\t"<< t3(Residual(i)) <<"\n";
  92.       cout << "\n\n";
  93.    }
  94.  
  95.    {
  96.       // traditional sum of squares and products method of calculation
  97.       // with subtraction of means - using Cholesky decomposition
  98.  
  99.       Matrix X(nobs,npred);
  100.       X.Col(1) << x1;  X.Col(2) << x2;
  101.       ColumnVector Y(nobs); Y << y;
  102.       ColumnVector Ones(nobs); Ones = 1.0;
  103.       RowVector M = Ones.t() * X / nobs;
  104.       Matrix XC(nobs,npred);
  105.       XC = X - Ones * M;
  106.       ColumnVector YC(nobs);
  107.       real m = Matrix(Ones.t() * Y)(1,1) / nobs;  YC = Y - Ones * m;
  108.       SymmetricMatrix SSQ; SSQ << XC.t() * XC;
  109.  
  110.       // Cholesky decomposition of SSQ
  111.       LowerTriangularMatrix L = Cholesky(SSQ);
  112.  
  113.       // calculate estimate
  114.       ColumnVector A = L.t().i() * (L.i() * (XC.t() * YC));
  115.  
  116.       // calculate estimate of constant term
  117.       real a = m - Matrix(M * A)(1,1);
  118.  
  119.       // Get variances of estimates from diagonal elements of invoice of SSQ
  120.       DiagonalMatrix D; D << L.t().i() * L.i();
  121.       ColumnVector V = D.CopyToColumn();
  122.       real v = 1.0/nobs + (L.i() * M.t()).SumSquare();
  123.  
  124.       // Calculate fitted values and residuals
  125.       ColumnVector Fitted = X * A + a;
  126.       ColumnVector Residual = Y - Fitted;
  127.       real ResVar = Residual.SumSquare() / (nobs-npred1);
  128.  
  129.       // print out answers
  130.       cout << "\nEstimates and their standard errors\n\n";
  131.       cout << a <<"\t"<< sqrt(v*ResVar) << "\n";
  132.       for (int i=1; i<=npred; i++)
  133.      cout << A(i) <<"\t"<< sqrt(V(i)*ResVar) << "\n";
  134.       cout << "\nObservations, fitted value, residual value\n";
  135.       for (i=1; i<=nobs; i++)
  136.      cout << X(i,1) <<"\t"<< X(i,2) <<"\t"<< Y(i) <<"\t"<<
  137.      t3(Fitted(i)) <<"\t"<< t3(Residual(i)) <<"\n";
  138.       cout << "\n\n";
  139.    }
  140.  
  141.    {
  142.       // Householder triangularisation method
  143.  
  144.       // load data - 1s into col 1 of matrix
  145.       Matrix X(nobs,npred1); ColumnVector Y(nobs);
  146.       X.Col(1) << 1.0;  X.Col(2) << x1;  X.Col(3) << x2;  Y << y;
  147.  
  148.       // do Householder triangularisation
  149.       // no need to deal with constant term separately
  150.       Matrix XT = X.t();             // Want data to be along rows
  151.       RowVector YT = Y.t();
  152.       LowerTriangularMatrix L; RowVector M;
  153.       HHDecompose(XT, L); HHDecompose(XT, YT, M); // YT now contains resids
  154.       ColumnVector A = L.t().i() * M.t();
  155.       ColumnVector Fitted = X * A;
  156.       real ResVar = YT.SumSquare() / (nobs-npred1);
  157.  
  158.       // get variances of estimates
  159.       L << L.i();
  160.       DiagonalMatrix D; D << L.t() * L;
  161.  
  162.       // print out answers
  163.       cout << "\nEstimates and their standard errors\n\n";
  164.       for (int i=1; i<=npred1; i++)
  165.      cout << A(i) <<"\t"<< sqrt(D(i)*ResVar) << "\n";
  166.       cout << "\nObservations, fitted value, residual value\n";
  167.       for (i=1; i<=nobs; i++)
  168.      cout << X(i,2) <<"\t"<< X(i,3) <<"\t"<< Y(i) <<"\t"<<
  169.      t3(Fitted(i)) <<"\t"<< t3(YT(i)) <<"\n";
  170.       cout << "\n\n";
  171.    }
  172. }
  173.  
  174. real t3(real r) { return int(r*1000) / 1000.0; }
  175.