regg.m

多维数据处理：MATLAB源程序用于处理多维数据

function [G,stdG]=regg(X,DimX,Fac,Factors,Weights);

ord=length(DimX);
if ord<3
   disp(' ')
   disp(' !!Corcondia only applicable for three- and higher-way arrays!!')
   break
end

if length(Fac)==1
   Fac=Fac*ones(1,ord);
end

vecX=X(:); % Vectorize X

% Make sure Weights are defined (as ones if none given)
if nargin<5
   Weights=ones(size(X));
end
if length(Weights(:))~=length(X(:));
   Weights=ones(size(X));
end
Weights=Weights(:);

% Set weights of missing elements to zero
id=find(isnan(vecX));
Weights(id)=zeros(size(id));
vecX(id)=zeros(size(id));

% Create Kronecker product of all but the last mode loadings
i1=sum(DimX(1:ord-2).*Fac(1:ord-2))+1;
i2=sum(DimX(1:ord-1).*Fac(1:ord-1));
L2=reshape(Factors(i1:i2),DimX(ord-1),Fac(ord-1));
i1=sum(DimX(1:ord-3).*Fac(1:ord-3))+1;
i2=sum(DimX(1:ord-2).*Fac(1:ord-2));
L1=reshape(Factors(i1:i2),DimX(ord-2),Fac(ord-2));
Z = kron(L2,L1);
for o=ord-4:-1:1
   i1=sum(DimX(1:ord-o-1).*Fac(1:ord-o-1))+1;
   i2=sum(DimX(1:ord-o).*Fac(1:ord-o));
   L=reshape(Factors(i1:i2),DimX(ord-o),Fac(ord-o));
   Z = kron(Z,L);
end

% Make last mode loadings, L
i1=sum(DimX(1:ord-1).*Fac(1:ord-1))+1;
i2=sum(DimX(1:ord).*Fac(1:ord));
L=reshape(Factors(i1:i2),DimX(ord),Fac(ord));

% We want to fit the model ||vecX - Y*vecG||, where Y = kron(L,Z), but 
% we calculate Y'Y and Y'vecX by summing over k
J=prod(DimX(1:ord-1));
Ytx = 0;
YtY = 0;
for k=1:DimX(ord)
   W=Weights((k-1)*J+1:k*J);
   WW=(W.^2*ones(1,prod(Fac)));
   Yk  = kron(L(k,:),Z);
   Ytx = Ytx + Yk'*(W.*vecX((k-1)*J+1:k*J));
   YtY = YtY + (Yk.*WW)'*Yk;
end

G=YtY\Ytx;

if nargout>1
   se = (sum(vecX.^2) + G'*YtY*G -G'*Ytx);
   mse = se/(length(vecX)-length(G));
   stdG=sqrt(diag(pinv(YtY))*mse);
end