parafacbig.m

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

% parafacbig.m.m
% uses cpfunc.m, etc. 
disp(' ');
disp(' CANDECOMP/PARAFAC program for three-way data with a big A-mode (I>>JK) ');
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('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(labb)~=m,disp(' Error: vector of labels for B (labb) has incorrect size');return;end;
if length(labc)~=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=[];

% To reduce computational cost, we follow the procedure by Harshman & Kiers (1997)
% to reduce the size for the array for the computational part

[Q,R]=qr(Xprep,0);
Xsmall=R;
n_original=n;
n=m*p;

% algorithm

if start==0
   [A,B,C,f,fp,iter,tripcos]=cpfunc(Xsmall,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(Xsmall,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;

% resize the A-mode 

A=Q*A;
n=n_original;

% ensure that column sums in A and B are positive

tau=diag(sign(sum(A)));
A=A*tau;
B=B*tau;
tau=diag(sign(sum(B)));
B=B*tau;
C=C*tau;

fprintf(' –andecomp/Parafac analysis with %g components gave a fit of %4.2f %% \n',r,fp); 
disp('  ');
disp(' The component matrices are normalized such that A and B have unit sums of squares');
disp(' To see component matrices, press "Enter" ');
pause;
disp(' ');
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 ');

imagedisp=input(' Do you want to display A-mode components as images? Press "1" ');
if isempty(imagedisp),imagedisp=0;end;
if imagedisp==1
   nx=input(' How many x-coordinates in each image: ');
   ny=input(' How many y-coordinates in each image: ');
   S=zeros(nx,ny);
   figure(1);
   for g=1:r,
      subplot(2,2,g)
      S(:)=A(:,g);
      imagesc(S);
   end;
end
   
disp('  ');