This routine will use Singular Value Decomposition and returns an expression sequence, thus two variables appear on the left hand side of the equals, one to hold the R matrix and the other the S matrix.
Helper Maple Code
Here are the various procedures that simplify the development of the examples below. This section is NOT intended to read unless one has an intimate interest in how the examples are built in Maple.
Library Calls and Global Parameters
This worksheet uses statistics, 3D plotting and the new Linear Algebra package
>
restart():
Digits := 32:
with(LinearAlgebra):
with(stats):
with(plots):
Warning, the name changecoords has been redefined
Here are some global parameters used to determine such things as the number of points to sample per class
>
pointsPerClass := 100:
animationSamples := 60:
Set Random Number Seed
It is handy to have the same point samples each time this file is run
> randomize(1111356902):
General Purpose Helper Procedures
Here are some procedures that are not directly related to Gaussian Random Variables, but which are helpful.
Reduce the precision of a Matrix to Aid Printing
This procedure simply rounds floating point numbers in Matrices to k digits so Maple pretty print looks better. This should not be necessary, but it appears to be.
>
sigDigits := proc(f,k)
local s, fn, tempDigits;
if (abs(f) < 10**(-k)) then
fn := 0.0;
else
s := ilog10(abs(f)) + 1;
# print("s is equal to ", s);
tempDigits := Digits;
Digits := max(s + k + 1,1);
# print("Setting Digits to ", Digits);
fn := (round (f*10**k))/(10.0**k);
Digits := tempDigits;
end if;
return fn;
end proc:
>
roundMatrix := proc(M,k)
map(proc(x) sigDigits(x,k) end proc,M);
end proc:
Generate and Manipulate a Data Matrix
This section includes procedures to generate sample data matrices.
A Data Matrix of Gaussian Random Variable Samples
The following set of procedures provide a general way of rotation, scaling and translating the points represented as columns in a homogeneous coordinates data matrix. Thus, the matrix is of size 1..4 by 1..pointsPerClass. The affine transformation construction procedures are general purpose. The support the final procedure, genGaussData, that takes 9 arguments. This procedure can construct point samples for any 3D Gaussian Random Variable by starting with samples from one with mean zero and standard deviation one. The 9 arguments specify the transformation from the zero-one reference frame to the final reference frame. They also relate directly to the paramters of the Gaussian pdf in the final reference frame. The first three are the standard deviation of a general form 3D Gaussian Random Variable when viewed in its principal coordinate representation. This is another way of saying these are the x, y and z standard deviation before it is rotated to a new orientation. The next three arguments are the translation away from the origin. These are just the x, y and z mean values for the pdf. The last three values specify the rotation of the distribution in Euler angles. In other words, rotation about the x, y and then z axis.
Generate a 4x4 Rotation Matrix
Given three Euler Angles, return a 4x4 homogenious rotation Matrix.
>
homRotation := proc(theta,phi,xi)
local Data, ROTx, ROTy, ROTz, R33, R44, i, j;
ROTx := Matrix([[1,0,0],[0,cos(theta),-sin(theta)],[0,sin(theta),cos(theta)]]):
ROTy := Matrix([[cos(phi),0,-sin(phi)],[0,1,0],[sin(phi),0,cos(phi)]]):
ROTz := Matrix([[cos(xi),-sin(xi),0],[sin(xi),cos(xi),0],[0,0,1]]):
R33 := map(eval,Multiply(Multiply(ROTz,ROTy),ROTx)):
R44 := Matrix(4,4,0.0): R44[4,4] := 1.0:
R44[1..3,1..3] := R33:
return R44:
end proc:
Generate a 4x4 Scale Matrix
Given an x, y and z axis scale factor, return a 4x4 homogenious scale Matrix
>
homScale := proc(sx, sy, sz)
local S:
S := Matrix(4,4,0.0):
S[1,1] := sx: S[2,2] := sy: S[3,3] := sz: S[4,4] := 1.0:
return S:
end proc:
Generate a 4x4 Translation Matrix
Given translations along the x, y an z axis, return a 4x4 homogenious translation Matrix
>
homTranslation := proc(tx, ty, tz)
local T:
T := Matrix(4,4):
T[1..4,1..4] := IdentityMatrix(4,4):
T[1,4] := tx:
T[2,4] := ty:
T[3,4] := tz:
return T:
end proc:
Generate a 4x4 Scale, Rotation, Translation Matrix
Given scale parameters, rotation parameters, and translation parameters, create a single 4x4 matrix that executes the scaling, followed by rotation, followed by translation.
>
homSRT := proc(sx,sy,sz,tx,ty,tz,theta,phi,xi)
local S, R, T:
S := homScale(sx,sy,sz):
R := homRotation(theta,phi,xi):
T := homTranslation(tx,ty,tz):
return Multiply(T,Multiply(R,S)):
end proc:
Generate a 4x4 homogenous Matrix from 3x3 specifying upper left
>
homPad := proc(M33)
local M44;
M44 := Matrix(4,4, IdentityMatrix(4,4));
M44[1..3,1..3] := M33;
return M44;
end proc:
Generate a Gaussiant Distributed data matrix with mean zero, standard deviation one
Using the global variable samplesPerClass, generate that many 3D points which are sampled from a simple 3D Gaussian where the mean is zero and standard deviation is one along each dimension
>
homGaussZeroOne := proc()
local n, D3, D4:
n := pointsPerClass:
D3 := Matrix(3,n, [[stats[random,normald[0.0,1.0]](n)],
[stats[random,normald[0.0,1.0]](n)],
[stats[random,normald[0.0,1.0]](n)]]);
D4 := Matrix(4,n,1):
D4[1..3,1..n] := D3:
return D4;
end proc:
Generate a scaled, rotated and then translated Gaussian data matrix
Starting with a data matrix of 3D Gaussian Random Variable sample with mean zero and standard deviation one for all three axes, scale the samples along each axis by sx, sy, sz, then rotate the samples by Euler angles theta, phi, xi, and finally translate the samples by tx, ty, tz.
>
genGaussData := proc(sx,sy,sz,tx,ty,tz,theta,phi,xi)
local D3, D4:
D4 := Multiply(homSRT(sx,sy,sz,tx,ty,tz,theta,phi,xi),homGaussZeroOne());
# D3 := SubMatrix(D4,[1..3],[1..ColumnDimension(D4)]);
return D4;
end proc:
Three classes from means, standard deviations and angles.
This procedure returns a three element list, where each element is the data matrix for a different class. The parameters for the class are specified in three matrices, MU, SD, and EA. These contain in each successive column, the means, standard deviations, and Euler angles for each for the three classes
>
listClasses := proc(MU, SD, EA)
local n, i, L:
n := ColumnDimension(MU);
L := []:
for i from 1 to n do
L := [op(L),genGaussData(SD[1,i],SD[2,i],SD[3,i],
MU[1,i],MU[2,i],MU[3,i],
EA[1,i],EA[2,i],EA[3,i])]:
od:
return L:
end proc:
Append Class Data Matrices to form one Common Data Matrix
This procedure will create a single common data matrix out of k distinct classes
>
columnAppendDataMatrices := proc(C)
local A, i, k, r, c, n:
k := nops(C);
n := ColumnDimension(C[1]):
r := RowDimension(C[1]):
A := Matrix(r,k*n,1.0):
for i from 1 to k do
c := (i-1)*n + 1;
A[1..r,c..c+n-1] := C[i]:
od:
return A:
end proc:
Sample Mean and Covariance Matrices
This section has two procedures, one to estimate the mean from a data matrix. The other to estimate the covariance matrix by computing the scatter matrix.
Compute the Mean Vector from a Data Matrix
This procedure sums the elements in the x, y and z dimensions of the data matrix and divides the totals by n.
>
dataMatrixMean := proc(A)
Vector(3,[seq(stats[describe,mean](convert(A,listlist)[i]),i=1..3)]);
end proc:
Compute the x, y and z ranges for a Data Matrix
Use the statistics package range facility to get the range of values along each axis. Return as a sequence of three range speciffications. Thus, calling order would list three range variables on the left side of the assignment statement.
>
dataMatrixRange := proc(A,dims)
seq(stats[describe,range](convert(A,listlist)[i]),i=1..dims);
end proc:
Sample Covariance from a Data Matrix
This procedure returns the 3x3 scatter matrix from a data matrix.
>
scatterMatrix := proc(A)
local n, MU, NT, M, CV, D3:
n := ColumnDimension(A):
MU := dataMatrixMean(A):
NT := homTranslation(-MU[1],-MU[2],-MU[3]):
D3 := SubMatrix(Multiply(NT,A),[1..3],[1..n]):
CV := Multiply(D3,Transpose(D3));
return CV;
end proc:
>
sampleCovariance := proc(A)
return scatterMatrix(A) . (1/(ColumnDimension(A) -1));
end proc:
Scatter Matrix within and Between Classes
>
scatterWBClass := proc(P)
local SW, SB, mu, mus, n;
SW := scatterMatrix(P[1]) + scatterMatrix(P[2]) + scatterMatrix(P[3]):
mus := [seq(dataMatrixMean(P[i]),i=1..3)]:
n := [seq(ColumnDimension(P[i]),i=1..3)];
mu := (1/eval(add(n[i],i=1..3))) . eval(add(n[i].mus[i],i=1..3)):
SB := eval(add(n[i].OuterProductMatrix(mus[i]-mu,mus[i]-mu),i=1..3)):
return SW, SB;
end proc:
Diagonilize a Covariance Matrix
This is a simple routine that conforms that yields a rotation and scale such that it represents the change in coordinates back to the PCA space for a given covariance (or scattter) matrix.
Return R and S Matrices given Sample Covariance/Scatter
Given the sample covariance or scatter matrix, give the diagonal decomposition such that
This routine will use Singular Value Decomposition and returns an expression sequence, thus two variables appear on the left hand side of the equals, one to hold the R matrix and the other the S matrix.
>
DiagonalizeRS := proc(omega)
local i, SV, RM, SM;
RM, SV := SingularValues(omega,output=['U','S'],outputoptions['U']=[datatype=float]):
SM := Matrix(RowDimension(omega),ColumnDimension(omega),0);
for i from 1 to RowDimension(omega) do SM[i,i] := sqrt(SV[i]); od;
return(RM,SM):
end proc:
Routines for Plotting Classes and Discriminants
Here are routines that will plot three data matrices of points along with 2 basis vectors. There are two routines, one that plots only the points and the other that plots the points and the discriminants. In both cases, the points are color coded by class and shown in a scatter plot. When shown, the discriminants are are drawn as line segments. The length of the basis vectors will be plus and minus one standard deviation. They are drawn with their origin at the centroid of the joint data matrix.
View Specification
Return a sequence of three ranges suitable for use as the view argument to the plotting routine such that x, y and z are centered on the mean with range sufficient to express largest variation.
>
jointDataView := proc(CC,dims)
local m, r, i, width, middle;
r := [dataMatrixRange(CC,dims)];
m := map(proc(x) op(2,x) - op(1,x) end proc,r);
width := max(seq(m[i],i=1..dims));
middle := map(proc(x) (op(1,x) + op(2,x))/2.0 end proc,r);
seq((middle[i]-(width/2))..(middle[i]+(width/2)),i=1..dims);
end proc:
Basis vector endpoints
Return two endpoints for each of the two Fisher Basis Vectors. As arguments, take in the rotation matrix, the scale matrix, and the joint data matrix. The joint data matrix is used to determine the joint mean. The resulting endpoints are encode in a 4x4 matrix, with the first two columns representing the start and end point of the first discriminant, the third and fourth the second. There are four rows because the endpoints are expressed in homogenious coordinates. This is in anticipation of our later desire to affine transform the space in which both discriminants and points are expressed.
>
basisVectorEndPoints := proc(R,S,P)
local delta1, delta2, mu, r33, result;
result := Matrix(4,4,1);
mu := dataMatrixMean(columnAppendDataMatrices(P));
# delta1 := Transpose(ScalarMultiply(R[1,1..3],S[1,1]));
# Changed scale to that of first vector for better viewing.
# delta2 := Transpose(ScalarMultiply(R[2,1..3],S[1,1]));
delta1 := ScalarMultiply(R[1..3,1],S[1,1]);
delta2 := ScalarMultiply(R[1..3,2],S[1,1]);
r33 := Matrix([(mu - delta1, mu + delta1, mu - delta2, mu + delta2)]);
result[1..3,1..4] := r33;
return result;
end proc:
Plot the points and basis vectors
Generate a scatter plot of the different classes expressed by the points in the data matrices contained in the list P.
>
pointsToSeq := proc(P,row)
local i, c, s;
c := ColumnDimension(P);
s := NULL;
for i from 1 to c do
s := s,[seq(P[j,i],j=1..row)];
od:
return s;
end proc:
>
pointPlotSpec := proc(P,row)
local i, lst,col,colors, cn;
colors := matrix(3,3,0.0);
colors[1,1] := 1.0; colors[2,2] := 1.0; colors[3,3] := 1.0;
lst := [];
for i from 1 to nops(P) do
col := COLOR(RGB,colors[1,i],colors[2,i],colors[3,i]);
cn := ColumnDimension(P[i]);
if (row = 2) then
lst := [op(lst),PLOT(POINTS(pointsToSeq(P[i],row)),SYMBOL(CIRCLE,10),col)];
else
lst := [op(lst),PLOT3D(POINTS(pointsToSeq(P[i],row)),SYMBOL(CIRCLE,10),col)];
fi;
od;
seq(lst[i],i=1..nops(P));
end proc:
>
basisPlotSpec := proc(BV)
local i, j, pts, lst, col, colors;
lst := [];
colors := matrix(3,2,1.0);
colors[1,1] := 0.0; colors[2,2] := 0.0;
j := 1;
for i from 1 to 3 by 2 do
pts := convert(Transpose(BV[1..3,i..(i+1)]),listlist);
col := COLOR(RGB,colors[1,j],colors[2,j],colors[3,j]);
j := j + 1;
lst := [op(lst),PLOT3D(CURVES(pts),col)];
od;
seq(lst[i],i=1..2);
end proc:
>
plotPoints := proc(P,theta,phi,axesType)
local PP, v;
PP := columnAppendDataMatrices(P);
v := [jointDataView(PP,3)];
display({pointPlotSpec(P,3)},insequence=false,axes=axesType,labels=[X,Y,Z],
axesfont=[TIMES,ROMAN,14],labelfont=[TIMES,ROMAN,16],
orientation=[theta,phi],scaling=CONSTRAINED,view=v);
end proc:
>
plotPoints2D := proc(P)
local PP, v;
PP := columnAppendDataMatrices(P);
v := [jointDataView(PP,2)];
display({pointPlotSpec(P,2)},insequence=false,axes=FRAMED,labels=[X,Y],
axesfont=[TIMES,ROMAN,14],labelfont=[TIMES,ROMAN,16],
scaling=CONSTRAINED,view=v);
end proc:
>
pointsToSeq := proc(P,row)
local i, c, s;
c := ColumnDimension(P);
s := NULL;
for i from 1 to c do
s := s,[seq(P[j,i],j=1..row)];
od:
return s;
end proc:
>
plotPointsBasis := proc(P,BV,theta,phi,axesType)
local PP, v, C1, C2, C3, B1, B2;
PP := columnAppendDataMatrices(P);
v := [jointDataView(PP,3)];
C1,C2,C3 := pointPlotSpec(P,3);
B1,B2 := basisPlotSpec(BV);
display({C1,C2,C3,B1,B2},insequence=false,axes=axesType,labels=[X,Y,Z],
axesfont=[TIMES,ROMAN,14],labelfont=[TIMES,ROMAN,16],
orientation=[phi,theta],scaling=CONSTRAINED,view=v);
end proc:
Compute R and S using the Generalized Eigenvector Method
>
generalizedEigenvectorsRS := proc(SB,SW)
local i, R, S, result:
result := Eigenvectors(SB,SW,output=list);
result := sort(result,proc(x,y) if (Re(x[1]) > Re(y[1])) then true else false fi end proc);
R := Matrix(3,3, 0.0);
S := Matrix(3,3, 0.0);
for i from 1 to 3 do
R[1..3,i] := result[i,3];
S[i,i] := sqrt(result[i,1]);
od;
return(map(Re,R),map(Re,S));
end proc:
Compute The Fisher Criterion Function
The criterion for the Fisher Basis Vectors maximizes a ration of two determinants. Here the arguement FB is a 3 row by 2 column matrix with one Fisher Basis Vector per column. SW is the within scatter matrix, and SB is the between class scatter matrix.
>
fischerCriterion := proc(FB,SW,SB)
local num, den, ret;
num := Multiply(Transpose(FB),Multiply(SB,FB));
den := Multiply(Transpose(FB),Multiply(SW,FB));
ret := Determinant(num) / Determinant(den);
convert(ret,float);
end proc:
Project onto 2D Fisher Basis Vectors
>
projectLD := proc(P,LD)
local Pp, LDP, PP;
LDP := homPad(LD)[1..4,1..2];
PP := [seq(Transpose(LDP).P[i],i=1..3)];
end proc: