Jonathan Bailleul schrieb:
>
>
> From a more global point of view, I would be interested in some "ready
> to use" source code/library to perform Procrustes analysis (I've heard
> about implementations named GPA). It seems like a reasonable scheme from
> the point of view of someone lacking the fundamental notions required to
> understand how it actually works...
>
> But of course I would be very pleased if you could "help me out" in
> explaining me, even very briefly, some of these notions.
>
Maybe, this code snippet is helpful? It's pascal/delphi.
AMAT is the matrix to be rotated columnwise; BMAT is the
destination matrix to what AMAT should be approximated.
The routine contains also a rowselector/colselector to
dynamically select only certain rows/columns for the computing
of the criteria and/or rotate only in selected columns.
Gottfried Helms
----------------------
procedure m_RotierenRahmen2(var amat:tmatrix;const bmat:tmatrix;const
rowselector,colselector:tintvector);
// the function-variable MyRotGet has to be set to an appropriate
// criteria-function before; i.e. PCA,Quartimax,Procrustes...
var
zl,i,s1,s2,s_1,s_2,spa:integer;
cs,sn,maxsn,mincs:real;
itercount,Maxiterate:integer;
begin
maxIterate:=500;
itercount:=0;
spa:=highsp(ColSelector);
for itercount:=1 to MaxIterate do
begin
maxsn:=0; mincs:=1;
for s_1:=1 to spa-1 do
for s_2:=s_1+1 to spa do
begin
s1:=Colselector[s_1]; // select a pair of columns for rotation
s2:=colSelector[s_2];
cs:=0;sn:=0;
if MyRotGet2(amat,bmat,s1,s2,cs,sn,RowSelector) // Getting the cosine&sine
for the rotation
// from loadings of all
(selected) rows
then begin
RotiereSpalten(amat,s1,s2,cs,-sn); // rotate this pair of
columns
if abs(sn)>maxsn then maxsn:=abs(sn);
if abs(cs)<mincs then mincs:=abs(cs);
end;
end;
if (maxsn<1e-16) and (abs(mincs-1)<1e-16) then break;
end;
end;
function RotKritGet_Another(var wmat:tmatrix;const bmat:tmatrix;s1,s2:integer;var
cs,sn:real;const RowSelector:tIntVector):boolean;
var
z,zl,z_:integer;
wx,wy,zx,zy,lae,c_su,s_su,c,s:real;
begin
result:=false;
c_su:=0; s_su:=0; // will hold the summed coefficents, giving the rotation angle
zl:=HighZl(Rowselector);
for z_:=1 to zl do begin z:=rowselector[z_]; // summing deviances over all rows
wx:=wmat[z,s1];
wy:=wmat[z,s2];
zx:=bmat[z,s1];
zy:=bmat[z,s2];
c:=wx*zx+wy*zy;
s:=wx*zy-wy*zx;
c_su:=c_su+c ;
s_su:=s_su+s;;
end;
cs:=c_su;sn:=s_su; // the summed coefficients over all(selected) rows
// represent the rotation-angle
lae:=sqrt(cs*cs+sn*sn);
if lae<1e-32 then exit;
cs:=cs/lae; sn:=sn/lae;
result:=true;
end;
.
.
=================================================================
Instructions for joining and leaving this list, remarks about the
problem of INAPPROPRIATE MESSAGES, and archives are available at:
. http://jse.stat.ncsu.edu/ .
=================================================================
