gemanova.m

进行方差分析的MATLAB源码

		for i=1:length(iteroff) % estimate each set of offsets one at a time
			%[OffsetsFromNstatReduced,OffsetsFromNstat] = nstat(xprep,'mean',iteroff{i},1);
			[OffsetsFromNstatReduced,OffsetsFromNstat] = nstat(xprep,'mean',iteroff{i});
			Offset{i}.OffsetsReduced = OffsetsFromNstatReduced;
			Offset{i}.Offsets = OffsetsFromNstat;
			xprep = xprep - OffsetsFromNstat;
		end
		% Replace unfolded data (xuf), so that the multilinear part is fitted to the residuals of the data subtracted the offsets
		if DoWeight & Iter > 1 % Then the data to fit is not the raw data in x but the modified data according to above in xuf
			xprep = reshape(xuf{1}',xsize);
		else
			xprep = x;
		end
		for i=1:length(iteroff) % estimate each set of offsets one at a time
			xprep = xprep - Offset{i}.Offsets;
		end
		for i = 1:order
			xuf{i} = unfoldmw(xprep,i)';
		end
	end
	
	for i = 1:order
		% Multiply the loads of all the orders together
		% except for the order to be estimated
		xuflo{i}=outerm(x0,i,1);
		
		% Regress the actual data on the estimate to get new loads in order i
		
		% Unconstrained
		if constraints(i) == 0
			ordiest = xuflo{i}\xuf{i};
			
			% nonnegativity   
		elseif constraints(i) == 1 
			ordiest = zeros(nocomp,xsize(i));
			xufloT = xuflo{i}'*xuflo{i}; % Calculate xproduct before loop
			for k = 1:xsize(i)
				ordiest(:,k) = CrossProdFastNnls(xufloT,xuflo{i}'*xuf{i}(:,k));
			end
			if any(sum(ordiest,2)==0);
				FeasibilityProblems=1;
				ordiest = .1*ordiest+.9*x0{i}';
			else
				FeasibilityProblems=0;
			end
			
			%Unimodality
		elseif constraints(i) == 2 
			ordiest=unimodal(xuflo{i},xuf{i},x0{i})';
			if any(sum(ordiest,2)==0);
				FeasibilityProblems=1;
				ordiest = .1*ordiest+.9*x0{i}';
			else
				FeasibilityProblems=0;
			end
			
			%Orthogonality
		elseif constraints(i) == 3
			% if this is the mode holding the scales modify so that loadings are not forced to be length one
			if i == ScaleMode
				Z = [];
				for fac = 1: nocomp
					Z = [Z kron(x0{i}(:,fac)/norm(x0{i}(:,fac)),xuflo{i}(:,fac))];
				end
				Scales = pinv(Z'*Z)*(Z'*xuf{i}(:));
			else
				Scales = ones(1,nocomp);
			end
			ZtX = (xuflo{i}*diag(Scales))'*xuf{i};
			ordiest=((ZtX*ZtX')^(-.5))*ZtX;
			if i == ScaleMode 
				ordiest = diag(Scales)*ordiest;
			end
			
			%Fixed
			
		elseif constraints(i) ~= -1 
			error([' The input constraints has not been correctly defined. The value ',num2str(constraints(i)),' is not possible'])
			
		elseif constraints(i) == -1 
			ordiest = x0{i}';
		end
		
		% Normalize the estimates (except the last (or other defined) order) and store them in the cell
		if i ~= ScaleMode & constraints(i)~=-1 % thus mode fixed
			Scal = 1./sqrt(sum(ordiest.^2,2));
			x0{i} = ordiest'*diag(Scal); % normalization leads to wrong model but that's corrected in the next update of the next mode, and for the last mode no normalization is performed, so that's ok, unless last mode is fixed.
			x0{ScaleMode} = x0{ScaleMode}*diag(Scal.^(-1));%ii = i+1;
		else
			x0{i} = ordiest';%ii = 1;
		end
		
	end
	% Calculate the estimate of the input array based on current loads
	xest = outerm(x0);
	% old version; xest = zeros(prod(xsize),nocomp);for j = 1:nocomp,      xvect = x0{1}(:,j);,      for ii = 2:order,         xvect = xvect*x0{ii}(:,j)';         xvect = xvect(:);      end,      xest(:,j) = xvect;,   end,   xest = sum(xest,2);   xest = reshape(xest,xsize);
	%Iterative preproc
	if iteroffsets
		XTrilin = xest; % Multilinear part to subtract from data before calculating offsets
		for ii=1:length(iteroff) % Add offsets one at a time
			xest = xest + Offset{ii}.Offsets;
		end      
	end
	xsq = xest;
	% Exchange missing with model estimates
	if Missing 
		x(MissId)=xest(MissId);
		for ii = 1:order
			xuf{ii} = unfoldmw(x,ii)';
		end
	end
	% Check to see if the fit has changed significantly
	if DoWeight
		xsq = ((x-xsq).*weights).^2;
	else
		xsq = (x-xsq).^2;
	end
	if Missing
		xsq(MissId)=0;  
	end
	ssq = sum(xsq(:));
	
	
	%disp(sprintf('On iteration %g ALS fit = %g',iter,ssq));
	if iter > 1 &~FeasibilityProblems
		abschange = abs(oldssq-ssq);
		relchange = abschange/ssq;
		if relchange < tol(1)
			flag = 1;
			if DumpToScreen
				disp(' '),disp('    Iteration    Rel. Change         Abs. Change         sum-sq residuals'),disp(' ')
				fprintf(' %9.0f       %12.10f        %12.10f        %12.10f    \n',iter,relchange,abschange,ssq);
				disp(' ')
				disp(' Iterations terminated based on relative change in fit error')
			end
		elseif abschange < tol(2)
			flag = 1;
			if DumpToScreen
				disp(' '),disp('    Iteration    Rel. Change         Abs. Change         sum-sq residuals'),disp(' ')
				fprintf(' %9.0f       %12.10f        %12.10f        %12.10f    \n',iter,relchange,abschange,ssq);
				disp(' ')
				disp(' Iterations terminated based on absolute change in fit error')
			end
		elseif iter > tol(3)-1
			flag = 1;
			if DumpToScreen
				disp(' '),disp('    Iteration    Rel. Change         Abs. Change         sum-sq residuals'),disp(' ')
				fprintf(' %9.0f       %12.10f        %12.10f        %12.10f    \n',iter,relchange,abschange,ssq);
				disp(' ')
				disp(' Iterations terminated based on maximum iterations')
			end
		end
	end
	if rem(iter,Show) == 0&DumpToScreen
		if iter == Show|rem(iter,Show*30) == 0
			disp(' '),disp('    Iteration    Rel. Change         Abs. Change         sum-sq residuals'),disp(' ')
		end
		fprintf(' %9.0f       %12.10f        %12.10f        %12.10f    \n',iter,relchange,abschange,ssq);
	end
	oldssq = ssq;
	
	% Every fifth iteration do a line search if ls == 1
	if (iter/5 == round(iter/5) & ls == 1)
		% Determine the search direction as the difference between the last two estimates
		for ij = 1:order
			searchdir{ij} = x0{ij} - oldx0{ij};
		end
		if iteroffsets
			for ij=1:length(iteroff)
				searchdirOffs{ij} = Offset{ij}.Offsets - OldOffset{ij}.Offsets;
			end
		end      
		% Initialize other variables required for line search
		testmod = x0;
		sflag = 0; 
		i = 0; 
		sd = zeros(10,1); 
		sd(1) = ssq;
		xl = zeros(10,1);
		while sflag == 0
			for k = 1:order
				testmod{k} = testmod{k} + (2^i)*searchdir{k};
			end
			if iteroffsets
				testoffsets = Offset;
				for j=1:length(Offset)
					testoffsets{j}.Offsets = testoffsets{j}.Offsets + (2^i)*searchdirOffs{j};
				end
			end
			% Calculate the fit error on the new test model
			xest = outerm(testmod);
			%Iterative preproc
			if iteroffsets
				XTrilin = xest; % TO SUBTRACT FROM DATA 
				for j=1:length(iteroff) % Add offsets one at a time
					xest = xest + testoffsets{j}.Offsets;
				end
			end
			if DoWeight
				xsq = ((x-outerm(testmod)).*weights).^2;
			else
				xsq = (x-outerm(testmod)).^2;
			end
			if Missing
				xsq(MissId)=0;  
			end
			% Save the difference and the distance along the search direction
			sd(i+2) = sum(xsq(:));
			xl(i+2) = xl(i+1) + 2^i;
			i = i+1;
			% Check to see if a minimum has been exceeded once two new points are calculated
			if i > 1 
				if sd(i+1) > sd(i)
					sflag = 1;
					% Estimate the minimum along the search direction
					xstar = sum((xl([i i+1 i-1]).^2 - xl([i+1 i-1 i]).^2).*sd(i-1:i+1));
					xstar = xstar/(2*sum((xl([i i+1 i-1]) - xl([i+1 i-1 i])).*sd(i-1:i+1)));
					% Save the old model and update the new one
					oldx0 = x0;
					for k = 1:order
						x0{k} = x0{k} + xstar*searchdir{k};
					end
					if iteroffsets
						OldOffset = Offset;
						for j=1:length(iteroff)
							Offset{j}.Offsets = Offset{j}.Offsets + + xstar*searchdirOffs{j};
						end
					end      
					
				end
			end
		end 
		% Calculate the estimate of the input array based on current loads
		xest = outerm(x0);
		%Iterative preproc
		if iteroffsets
			XTrilin = xest; % TO SUBTRACT FROM DATA 
			for j=1:length(iteroff) % Add offsets one at a time
				xest = xest + Offset{j}.Offsets;
			end      
		end
		xsq = xest;
		if DoWeight
			xsq = ((x-xsq).*weights).^2;
		else
			xsq = (x-xsq).^2;
		end
		if Missing
			xsq(MissId)=0;  
		end
		oldssq = sum(xsq(:));
		if Missing 
			x(MissId)=xest(MissId);
			for j = 1:order
				xuf{j} = unfoldmw(x,j)';
			end
		end
		% disp(sprintf('SSQ at xstar is %g',oldssq))
	else
		% Save the last estimates of the loads
		oldx0 = x0;
	end
	% Exchange missing with model estimates
end


% Plot the loadings
if plots ~= 0
	h1 = figure('position',[170 130 512 384],'name','PARAFAC Loadings');
	for i = 1:order
		subplot(order,1,i)
		if ~isempty(scl{i})
			if xsize(i) <= 50
				plot(scl{i},x0{i},'-+')
			else
				plot(scl{i},x0{i},'-')
			end
		else
			if xsize(i) <= 50
				plot(x0{i},'-+')
			else
				plot(x0{i},'-')
			end
		end
		ylabel(sprintf('Dimension %g',i))
		if i == 1
			title('Loadings for Each Dimension')
		end
	end
end

% Calculate and plot the residuals  
dif = (x-outerm(x0)).^2;
res = cell(1,order);
for i = 1:order
	x = dif;
	for j = 1:order
		if i ~= j
			x = sum(x,j);
		end
	end
	x = squeeze(x);
	res{i} = x(:);
	if plots ~= 0
		if i == 1  
			figure('position',[145 166 512 384],'name','PARAFAC Residuals')
		end
		subplot(order,1,i)
		if ~isempty(scl{i})
			if xsize(i) <= 50
				plot(scl{i},res{i},'-+')
			else
				plot(scl{i},res{i},'-')
			end
		else
			if xsize(i) <= 50
				plot(res{i},'-+')
			else
				plot(res{i},'-')
			end
		end
		ylabel(sprintf('Dimension %g',i))
		if i == 1
			title('Residuals for Each Dimension')
		end
	end
end

% Bring the loads back to the front
if plots ~= 0
	figure(h1)
end

% Save the model as a structured array    
model = struct('xname',inputname(1),'name','PARAFAC','date',date,'time',clock,...
	'size',xsize,'nocomp',nocomp,'tol',tol,'final',[relchange abschange iter]);
model.ssq = [tssq ssq];
model.loads = x0;
model.res = res;
model.scl = scl;
%Iterative preproc
if iteroffsets
	for i=1:length(iteroff) % Add offsets one at a time
		model.Offsets.Parameters{i} = Offset{i}.OffsetsReduced;
		model.Offsets.Model{i} = Offset{i}.Offsets;
		model.Offsets.Def{i} = iteroff{i};
	end      
end

function B=unimodal(X,Y,Bold)

% Solves the problem min|Y-XB'| subject to the columns of 
% B are unimodal and nonnegative. The algorithm is iterative
% If an estimate of B (Bold) is given only one iteration is given, hence
% the solution is only improving not least squares
% If Bold is not given the least squares solution is estimated
%
% Copyright 1997
%
% Rasmus Bro
% Royal Veterinary & Agricultural University
% Denmark
% [e-mail]
%
% Reference
% Bro and Sidiropoulos, "Journal of Chemometrics", 1998, 12, 223-247. 


if nargin==3
	B=Bold;
	F=size(B,2);
	for f=1:F
		y=Y-X(:,[1:f-1 f+1:F])*B(:,[1:f-1 f+1:F])';
		beta=pinv(X(:,f))*y;
		B(:,f)=ulsr(beta',1);
	end
else
	F=size(X,2);
	maxit=100;
	B=randn(size(Y,2),F);
	Bold=2*B;
	it=0;
	while norm(Bold-B)/norm(B)>1e-5&it<maxit
		Bold=B;
		it=it+1;
		for f=1:F
			y=Y-X(:,[1:f-1 f+1:F])*B(:,[1:f-1 f+1:F])';
			beta=pinv(X(:,f))*y;
			B(:,f)=ulsr(beta',1);
		end
	end
	if it==maxit
		disp([' UNIMODAL did not converge in ',num2str(maxit),' iterations']);
	end
end


function [b,All,MaxML]=ulsr(x,NonNeg);

% ------INPUT------
%
% x          is the vector to be approximated
% NonNeg     If NonNeg is one, nonnegativity is imposed
%
%
%
% ------OUTPUT-----
%
% b 	     is the best ULSR vector
% All 	     is containing in its i'th column the ULSRFIX solution for mode
% 	     location at the i'th element. The ULSR solution given in All
%            is found disregarding the i'th element and hence NOT optimal
% MaxML      is the optimal (leftmost) mode location (i.e. position of maximum)
%
% ___________________________________________________________
%
%
%               Copyright 1997
%
% Nikos Sidiroupolos
% University of Maryland
% Maryland, US
%
%       &
%
% Rasmus Bro
% Royal Veterinary & Agricultural University
% Denmark
%
% 
% ___________________________________________________________


% This file uses MONREG.M

x=x(:);
I=length(x);
xmin=min(x);
if xmin<0
	x=x-xmin;
end


% THE SUBSEQUENT 
% CALCULATES BEST BY TWO MONOTONIC REGRESSIONS

% B1(1:i,i) contains the monontonic increasing regr. on x(1:i)
[b1,out,B1]=monreg(x);

% BI is the opposite of B1. Hence BI(i:I,i) holds the monotonic
% decreasing regression on x(i:I)
[bI,out,BI]=monreg(flipud(x));
BI=flipud(fliplr(BI));