Category : C++ Source Code
Archive   : YAMP15.ZIP
Filename : APP2.TEX

\chapter{testreg.cpp}
This is the test program for the regression decompositions
in the VMatrix system. The singular valued decomposition
inverts the 11x11 Hilbert matrix with a maximum error of 0.001.

\begin{verbatim}


#include "virt.h"

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

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

VMatrix& ExactHilbInv( int n )
{
// construct exact hilbert matrix inverse
Dispatch->Inclevel();
VMatrix Hi("Hilbert Inverse",n,n);

double diag, dn = n;
for( int i=1; i<=n; i++){
double im1 = (double) (i-1.0);
diag = ( i == 1 ) ? dn : ((dn-im1)*diag*(dn+im1))/(im1*im1);

long double rest = diag*diag;
Hi.M(i,i) = rest/(2.0*im1+1.0);
for( int j=i+1; j<=n; j++ ){
double jm1 = (double) (j-1);
rest = -((dn-jm1)*rest*(dn+jm1))/(jm1*jm1);
Hi.M(i,j) = rest/(im1+jm1+1.0);
Hi.M(j,i) = Hi.m(i,j);
}
}
Dispatch->Push( Hi );
return Dispatch->DecReturn();
}


void TestHilb( int hilb_ind )
// test svd on hilbert matrix
{
Dispatch->Inclevel();
VMatrix hilb = Kron( Fill(1,hilb_ind,1), Index(1, hilb_ind));
// h(i,j) = 1.0/(i+j-1)
hilb = 1.0 / ( (hilb+Tran(hilb))-1.0);
hilb.Nameit("Hilbert Matrix");
hilb.DisplayMat();
// use formula. Looses accuracy of 0.001 at n>=11
(ExactHilbInv( hilb_ind )*hilb ).DisplayMat();

// use sweep
VMatrix Gauss = Inv(hilb)*hilb;
Gauss.Nameit("Gaussian elimination");
Gauss.DisplayMat();

// use singular valued decomposition
VMatrix S, V, D;
Svd( hilb, S, V, D);
VMatrix 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) = 1.0 / tt;
}

VMatrix Ggauss = (D*t*Tran(S))*hilb;
Ggauss.Nameit("Generalized inversion");
Ggauss.DisplayMat();

Dispatch->Declevel();
}

///////////////////// now test regressions on sample data

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

VMatrix x = Index( 1, N ) - N/2;
VMatrix c1 = Fill(N, 1, 1.0);
VMatrix 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 = 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 ®ression(VMatrix& x, VMatrix& y)
// do a multiple linear regression
{
Dispatch->Inclevel();

int N = x.r, p = x.c;
// solve for regression parameters using sweep
VMatrix yx = Ch(y, x);
VMatrix 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
VMatrix 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;

QR(x, Q, R);
VMatrix 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;

Svd(x, S, V, D);
VMatrix 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) <= 0.0) ? 0.0 : 1.0 / tt;
}
VMatrix betahat = D*t*Tran(S)*y;

Dispatch->Push(betahat);
return Dispatch->DecReturn();
}
VMatrix &GetSerialCovar(VMatrix &R)
{
Dispatch->Inclevel();

double n = (double) R.r;
VMatrix centered = R - Sum(R / n).m(1, 1);
VMatrix z = Fft(centered);
VMatrix spectdens = Sum(z % z / n, COLUMNS);
VMatrix covar = Fft(spectdens, FFALSE);

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

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

TestHilb( 11 );

Setwid(15);
Setdec(10);


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

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();

Setwid(6);
Setdec(3);

vclose();
}

\end{verbatim}