\chapter{testreg.cpp}
This is the test program for the regression decompositions 
in the VMatrix system.

\begin{verbatim}

#include "virt.h"

//required global declaration for the
//  matrix stack object

unsigned int _stklen = STACKLENGTH;
MStack *Dispatch = new MStack;


VMatrix &getx( int N )
  // create an x matrix
  {
      Dispatch->Inclevel();
      VMatrix x, c1, x2;

      x = Fill(N,1,0);
      c1 = Fill(N,1,1.0);
      for( int i=1; i<=N; i++)
         x.M(i,1) = (double) (i-10);
      x2 = x%x;
      x = Ch( c1, Ch( x, x2) );

      // push x onto the stack
      Dispatch->Push(x);
      // decrement the subroutine nesting level
      // and return the stack top
      return Dispatch->DecReturn();
  }

VMatrix &gety( VMatrix &x, VMatrix &beta)
  // create a y matrix
  {
      Dispatch->Inclevel();
      VMatrix y;

      y = x*beta;
      srand(123);
      for(int i=1; i<=y.r; i++) {
        // use sum of 3 uniforms for an approximate
        // normal random variable
        int u = random(100)+random(100)+random(100)+3;
        y.M(i,1) = y.m(i,1) + ((double) (u-150))/300.0;
      }
      Dispatch->Push(y);
      return Dispatch->DecReturn();
  }

VMatrix &regression( VMatrix& x, VMatrix& y )
// do a multiple linear regression
      {
      Dispatch->Inclevel();
      VMatrix yx, reg, betahat;
      int N=x.r, p=x.c;

      // solve for regression parameters using sweep
      yx = Ch(y,x);
      reg = Sweep( 2,p+1, Tran(yx)*yx);
      // calculate mean square error
      reg.M(1,1) = reg.m(1,1)/((double )( N-p));
      reg.DisplayMat();


      // solve regression using normal equations
      betahat = Inv(Tran(x)*x)*Tran(x)*y;

      Dispatch->Push( betahat );
      return Dispatch->DecReturn();
}
VMatrix &QRregression( VMatrix &x, VMatrix &y)
  {
     // use QR decomposition to solve regression

     Dispatch->Inclevel();
     VMatrix Q, R, betahat;

     QR( x, Q, R);
     betahat = Inv( Tran(R)*R )*Tran(R)*Tran(Q)*y;

     Dispatch->Push( betahat );
     return Dispatch->DecReturn();
  }
VMatrix &GinvRegression( VMatrix &x, VMatrix &y )
  {
     // use Ginv to solve normal equations
     Dispatch->Inclevel();

     VMatrix betahat = Ginv( Tran(x)*x )*Tran(x)*y;

     Dispatch->Push( betahat );
     return Dispatch->DecReturn();
  }


VMatrix &SVDregression( VMatrix &x, VMatrix &y)
  {
     // use SVD to solve regression x = S Diag(V) Tran(D)
     Dispatch->Inclevel();
     VMatrix S, V, D, betahat, t;


     Svd( x, S, V, D);
     t = Fill( V.r, V.r, 0.0);
     for( int i=1; i<=V.r; i++ ){
       double tt = V.m(i,1);
       t.M(i,i) = (fabs(tt) <= 1.0e-10 ) ? 0.0 : 1.0/tt;
     }
     betahat = D*t*Tran(S)*y;
     Dispatch->Push( betahat );
     return Dispatch->DecReturn();
  }
VMatrix &GetSerialCovar( VMatrix &R )
  {
    Dispatch->Inclevel();
    VMatrix centered, z, spectdens, covar;
    double n = (double) R.r;
    centered = R - Sum( R/n ).m(1,1);
    z = Fft( centered );
    spectdens = Sum( z%z/n, COLUMNS);
    covar = Fft( spectdens, FFALSE );
    Dispatch->Push( covar );
    return Dispatch->DecReturn();
  }

main()
{
      Dispatch->Inclevel();
      VMatrix x, beta("beta",3,1), y, betahat, resids, serial;
      Setwid(15);
      Setdec(10);


      beta.M(1,1) = 1;
      beta.M(2,1) = 0.5;
      beta.M(3,1) = 0.25;

      // Also try this with 200, and note how poorly the
      // ginv regression does because of colinearities

      x = getx( 55 );
      y = gety(x,beta);

      betahat = regression(x,y);
      betahat.Nameit( "Text book betahat");
      (Tran(beta)).DisplayMat();
      (Tran(betahat)).DisplayMat();

      betahat = QRregression( x,y );
      betahat.Nameit( "QR betahat");
      (Tran(beta)).DisplayMat();
      (Tran(betahat)).DisplayMat();

      betahat = GinvRegression( x, y);
      betahat.Nameit( "Ginv betahat");
      (Tran(beta)).DisplayMat();
      (Tran(betahat)).DisplayMat();

      betahat = SVDregression( x,y);
      betahat.Nameit("SVD regression");
      (Tran(beta)).DisplayMat();
      (Tran(betahat)).DisplayMat();

      resids = y - x*betahat;
      serial = GetSerialCovar( resids );
      (Tran( Submat( serial, 5,1 ) )).DisplayMat();

   vclose();
}

\end{verbatim}