parafac.m

这是一个三面阵列的数学工具软件

% parafac.m
% uses cpfunc.m, etc. 
disp(' ');
disp(' Simple CANDECOMP/PARAFAC program');
disp(' ');
disp(' Henk A.L. Kiers, April 2001');
disp(' ');
disp(' Data expected on X; rest interactive');
disp(' ');
cc=0;
while cc==0
  n=input(' Specify the number of A-mode entities ');
  m=input(' Specify the number of B-mode entities ');
  p=input(' Specify the number of C-mode entities ');
  if isempty(n)|isempty(m)|isempty(p),cc=0;disp(' Please specify ALL sizes! ');
  else,cc=1;end;
end;

if exist('X')==0,disp(' Error: Data not available on X');return;end;
Xprep=X;
[rr,cc]=size(Xprep);
if rr~=n | cc~=m*p, disp(' Error: array size does not agree with input data matrix X');
   return;
end;
if exist('laba')==0,laba=num2str([1:n]');aa=length(laba(1,:));laba(:,2:aa+1)=laba;laba(:,1)='A';end
if exist('labb')==0,labb=num2str([1:m]');aa=length(labb(1,:));labb(:,2:aa+1)=labb;labb(:,1)='B';end
if exist('labc')==0,labc=num2str([1:p]');aa=length(labc(1,:));labc(:,2:aa+1)=labc;labc(:,1)='C';end
if length(laba(:,1))~=n,disp(' Error: vector of labels for A (laba) has incorrect size');return;end;
if length(labb(:,1))~=m,disp(' Error: vector of labels for B (labb) has incorrect size');return;end;
if length(labc(:,1))~=p,disp(' Error: vector of labels for C (labc) has incorrect size');return;end;

disp('  ');
% Centering
Xprep=X;
disp('  ');
disp(' How do you want to center your array ?');
disp(' 1= across A-mode ');
disp(' 2= across B-mode ');
disp(' 3= across C-mode ');
disp(' 12= across A-mode and across B-mode ');
disp(' 13= across A-mode and across C-mode ');
disp(' 23= across B-mode and across C-mode ');
cc=input(' Specify your choice:  ');
if isempty(cc);cc=0;end;
centopt=cc;
disp('  ');
if cc==1 | cc==12 | cc==13
   Xprep=cent3(Xprep,n,m,p,1);
   disp(' Data have been centered across A-mode');
end;
if cc==2 | cc==12 | cc==23
   Xprep=cent3(Xprep,n,m,p,2);
   disp(' Data have been centered across B-mode');
end;
if cc==3 | cc==13 | cc==23
   Xprep=cent3(Xprep,n,m,p,3);
   disp(' Data have been centered across C-mode');
end;

% Normalizing 
disp('  ');
disp(' How do you want to normalize your array ?');
disp(' 1= within A-mode ');
disp(' 2= within B-mode ');
disp(' 3= within C-mode ');
cc=input(' Specify your choice:  ');
if isempty(cc);cc=0;end;
normopt=cc;
disp('  ');
if cc==1 
   Xprep=norm3(Xprep,n,m,p,1);
   disp(' Data have been normalized within A-mode');
end;
if cc==2 
   Xprep=norm3(Xprep,n,m,p,2);
   disp(' Data have been normalized within B-mode');
end;
if cc==3
   Xprep=norm3(Xprep,n,m,p,3);
   disp(' Data have been normalized within C-mode');
end;

disp('  ');
disp(' NOTE: ');
disp('   The preprocessed data are now available in Xprep');
disp('   Xprep can be used for analyses outside the Parafac.m program');
disp('  ');
disp(' ');
r=input(' How many components do you want to use?  ');
disp(' ');
conv=input(' Specify convergence criterion (default=1e-6)');
if isempty(conv),conv=1e-6;end;
disp(' ');
maxit=input(' Specify the maximum number of iterations you allow (default=5000)?  ');
if isempty(maxit),maxit=5000;end;
disp(' ');
cc=0;
cc=input(' Do you want to use constraints? If so, enter "1": ');
if isempty(cc),cc=0;end;
if cc==0
   ort1=1;ort2=1;ort3=1;disp(' No constraints imposed');
end;
if cc==1   
  disp(' 1= No constraints (default)');
  disp(' 2= Orthogonality constraints');
  disp(' 3= Zero correlations constraints');
  ort1=input(' Specify the A-mode constraints: ');
  ort2=input(' Specify the B-mode constraints: ');
  ort3=input(' Specify the C-mode constraints: ');
  if isempty(ort1),ort1=1;end;
  if isempty(ort2),ort2=1;end;
  if isempty(ort3),ort3=1;end;
  if ort1~=1&ort1~=2&ort1~=3,ort1=1;end;  
  if ort2~=1&ort2~=2&ort2~=3,ort2=1;end;  
  if ort3~=1&ort3~=2&ort3~=3,ort3=1;end;  
  if ort1==2 & r>n,
     fprintf('\n Warning : due to the constraint on A, you cannot use more than %g components \n',n);
     r=n;
  end;
  if ort2==2 & r>m,
     fprintf('\n Warning: due to the constraint on B, you cannot use more than %g components \n',m);
     r=m;
  end;
  if ort3==2 & r>p,
     fprintf('\n Warning: due to the constraint on C, you cannot use more than %g components \n',p);
     r=p;
  end;
  if ort1==3 & r>n-1,
     fprintf('\n Warning: due to the constraint on A, you cannot use more than %g components \n',n-1);
     r=n-1;
  end;
  if ort2==3 & r>m-1,
     fprintf('\n Warning: due to the constraint on B, you cannot use more than %g components \n',m-1);
     r=m-1;
  end;
  if ort3==3 & r>p-1,
     fprintf('\n Warning: due to the constraint on C, you cannot use more than %g components \n',p-1);
     r=p-1;
  end;
end;
disp('  ');
disp(' By default, only a rationally started analysis run will be carried out.');
disp(' To decrease the chance of missing the optimal solution, you may use');
disp(' additional, randomly started runs.');
addanal=input(' If you want additional runs, specify how many (e.g., 4):  ');
if isempty(addanal),addanal=0;end;
disp('  ');
start=0;
A=[];B=[];C=[];
if start==0
   [A,B,C,f,fp,iter,tripcos]=cpfunc(Xprep,n,m,p,r,ort1,ort2,ort3,start,conv,A,B,C,maxit);
end;
func=ones(1+addanal,1);
func=fp;
for run=1:addanal
   disp('  ');
   fprintf(' Run no. %g \n',run+1);
   disp('  ');
   start=1;
   [Ar,Br,Cr,fr,fpr,iterr,tripcosr]=cpfunc(Xprep,n,m,p,r,ort1,ort2,ort3,start,conv,A,B,C,maxit);
   func(run+1)=fpr;
   if fpr>1.0001*fp % if fit more than .01% better is found, replace solution
      A=Ar;B=Br;C=Cr;f=fr;iter=iterr;fp=fpr;tripcos=tripcosr;
   end;
   disp('  ');
end;
disp('  ');
if addanal>=1
   disp('  ');
   disp(' Fit % values from all runs:');
   disp(' ');
   writescr(func','6.2');
   disp('  ');
end;
% normalize solution
ssa=diag(sum(A.^2).^(-.5));
ssb=diag(sum(B.^2).^(-.5));
A=A*ssa;
B=B*ssb;
C=C/ssa/ssb;

% manual permutation and reflection:
disp('  ');
disp('  You can now manually PERMUTE and REFLECT columns/rows of solution');;
disp('  ');
cc=input(' If you want to reflect/permute columns/rows, specify "1":  ');
if isempty(cc);cc=0;end;
while cc==1
   disp('Here you find (max. first 10) rows of A, B and C');
   kk=min(n,10); disp(' A  ');writescr(A(1:kk,:),'6.2');
   kk=min(m,10); disp(' B  ');writescr(B(1:kk,:),'6.2');
   kk=min(p,10); disp(' C  ');writescr(C(1:kk,:),'6.2');
   tau=[];pp=[];
   disp(' ');
   tau=input(' Give vector for reflection of columns of A (e.g., [1 -1 -1 1 ..]) ');
   if isempty(tau),tau=ones(1,r);end;
   if ssq(tau)~=r,disp(' incorrect reflection');tau=ones(1,r);end;
   if ssq(tau.^2-1)~=0,disp(' incorrect reflection');tau=ones(1,r);end;
   disp(' ');
   pp=input(' Give vector with new order of columns of A, B and C  (e.g., [3 1 4 2 ..]) ');
   if isempty(pp),pp=1:r;end;
   if ssq(size(pp)-size(ones(1,r)))>0,disp(' incorrect permutation');pp=1:r;end;
   if ssq(sort(pp')'-[1:r])>0,disp(' incorrect permutation');pp=1:r;end;
   A=A*diag(tau);B=B*diag(tau);
   A=A(:,pp);B=B(:,pp);C=C(:,pp);
   disp('  ');disp(' B  ');disp('  ');writescr(B,'6.2')
   disp('  ');disp(' C  ');disp('  ');writescr(C,'6.2')
   tau=[];pp=[];
   disp(' ');
   tau=input(' Give vector for reflection of columns of B (e.g., [1 -1 -1 1 ..]) ');
   if isempty(tau),tau=ones(1,r);end;
   if ssq(tau)~=r,disp(' incorrect reflection');tau=ones(1,r);end;
   if ssq(tau.^2-1)~=0,disp(' incorrect reflection');tau=ones(1,r);end;
   disp(' ');
   B=B*diag(tau);C=C*diag(tau);
   disp('  ');
   disp('  ');disp(' C ');disp('  ');writescr(C,'6.2')
   tau=[];pp=[];
   disp(' ');
   tau=input(' Give vector for reflection of columns of C (e.g., [1 -1 -1 1 ..]) ');
   if isempty(tau),tau=ones(1,r);end;
   if ssq(tau)~=r,disp(' incorrect reflection');tau=ones(1,r);end;
   if ssq(tau.^2-1)~=0,disp(' incorrect reflection');tau=ones(1,r);end;
   disp(' ');
   C=C*diag(tau);A=A*diag(tau);
   disp('  ');disp(' C  ');disp('  ');writescr(C,'6.2')
   disp(' ');
   disp(' ');
   cc=input(' If you want to further reflect/permute columns/rows, specify "1":  ');
   if isempty(cc);cc=0;end;
end;


fprintf(' –andecomp/Parafac analysis with %g components gave a fit of %4.2f %% \n',r,fp); 
disp('  ');
disp(' Simple check on degeneracy: inspect matrix of triple congruences');
disp('  ');
writescr(phi(A,A).*phi(B,B).*phi(C,C),'8.4');
disp('  ');

disp(' The component matrices are normalized such that A and B have unit sums of squares');
disp(' To see component matrices (A only if n_A<40), press "Enter" ');
pause;
disp(' ');
if n<40,disp('  ');disp('  A  ');disp('  ');writescr(A,'6.2',laba);end
disp('  ');disp('  B  ');disp('  ');writescr(B,'6.2',labb);
disp('  ');disp('  C  ');disp('  ');writescr(C,'6.2',labc);
disp('  ');
disp(' The results are available in A, B, and C ');