gemanova.m

进行方差分析的MATLAB源码

 
 
 
 
 
 
 function model = parafac(x,nocomp,scl,tol,x0,options,plots,constraints,weights);

%PARAFAC Parallel factor analysis for n-way arrays
%  PARAFAC will decompose an array of order n (where n >= 3)
%  into the summation over the outer product of n vectors.
%  Missing values must be NaN or Inf.
%  
%  INPUTS
%  x           : the multi-way array to be decomposed
%  nocomp      : the number of components to estimate, 
%
%  OPTIONAL INPUTS
%  scl         : optional inputs of a cell array of vectors for plotting
%                loads against (set scl=[] for none, if elements of scl are
%                empty or of incorrect length they are reset),
%  tol         : convergence tolerance (a 1 by 3 vector consisting of 
%                the relative change in fit, absolute change in fit 
%                and maximum iterations {default tol = [1e-6 1e-6 10000]}), 
%  x0          : initial estimate of the factors, default is none = [].
%                x0 is a cell array of the same form as the fitted loadings
%                output in model.loads. Thus, x0{i} holds the loading matrix
%                in the i'th mode.
%                   If x0 is correctly provided, it will be used for initializing
%                the algorithm. Note, that if constraints(i) is set to -1, then
%                the loadings of that mode will not be updated, hence the initial
%                stay fixed. This way, one may use the algorithm for fitting an 
%                existing model to new data, by inputting the old model (model.loads) 
%                and fixing all except the sample mode loadings.
%                   It is also possible to input a complete PARAFAC model (a 
%                structure). This is useful when a prior model also includes
%                offsets, so that these are also used as initial values.
%  options     : A vector defining algorithmic settings:
%                options(1) = a flag which turns off the line search options when 0
%                options(2) = turns off screen output when set to zero
%                Default, options = [] = [1 1];
%  plots       : a flag which turns off the plotting of the loads and 
%                residuals when set to zero. 
%  constraints : Vector defining constraints. constraints has length equal to the number
%                of modes (for a three-way array, constraints is a 3-vector). 
%                constraints(1) =  0 => First mode loadings unconstrained
%                constraints(1) =  1 => First mode loadings nonnegativity-constrained
%                constraints(1) =  2 => First mode loadings unimodality-unconstrained (and nonneg)
%                constraints(1) =  3 => First mode loadings orthogonality-unconstrained
%                constraints(1) = -1 => First mode loadings fixed (requires initial values given)
%                constraints(2) =  ...  Defines constraints on second mode loadings
%                Example.: constraints = [0 1 2] means A unconstrained, B nonneg, & C unimodal
%  weights     : Enables weighted least squares fitting. Type <parafac('weights')> to get help on 
%                this feature
%
%  The output is a structured array (model) containing the 
%  PARAFAC model elements:
%
%     xname: name of the original workspace input variable
%      name: type of model, always 'PARAFAC'
%      date: model creation date stamp
%      time: model creation time stamp
%      size: size of the original input array
%     ncomp: number of components estimated
%       tol: tolerance vector used during model creation
%     final: final tolerances at termination
%       ssq: total and residual sum of squares
%     loads: 1 by order cell array of the loadings in each dimension
%       res: 1 by order cell array squared residuals summed over each dimension
%       scl: 1 by order cell array with scales for plotting loads 
%
%  This routine uses alternating least squares (ALS) in combination with
%  a line search every fifth iteration. For 3-way data, the intial estimate
%  of the loadings is obtained from the tri-linear decomposition (TLD).
%
%I/O: model = parafac(x,nocomp,scl,tol,x0,options,plots,constraints,weights,iteroff);
%
%See also:  GRAM, MPCA, MWFIT, OUTER, OUTERM, TLD, UNFOLDM, XPLDST

%Copyright Eigenvector Research, Inc. 1998-2000
%bmw March 4, 1998
%rb September, 1999 (added nnls, ulsr, fixed modes, orthogonality, weights, missing data, screendumps, moving the mode holding the scale in case the last mode is fixed)
%nbg 3/1/00 (check scl so routine doesn't bomb when not correct)
%rb March , 2000 (handle empty input for weights and constraints, do iteratively computed offsets)
%Ditto (included gui for help on iterative inputs)
%Ditto (corrected error so that old loading-values are accepted as cell arrays)
%rb April, 2000, corrects so that offsets can be included in old-values input, x0, by allowing a model struct to be input
%rb April, 2000, Fixed updating of unfolded X in linesearch when missing
%rb August, 2000, Added ATLD for initialization
%rb December, 2000, Fixed error in fixed mode loadings estimation
%rb February, 2001, removed iteroff (iterative offsets)


if nargin == 0
	disp(' ')
	disp(' PARAFAC Parallel factor analysis for n-way arrays - short help')
	disp(' ')
	disp(' I/O: model = parafac(x,nocomp,scl,tol,x0,options,plots,constraints,weights,iteroff);')
	disp(' ')
	disp(' Type <<help parafac>> for extended help')
	disp(' ')
	break
end

if nargin ==1
	if strcmp(lower(x),'weights')
		disp(' USING A WEIGHTED LOSS FUNCTION')
		disp(' ')
		disp(' Through the use of the input ''weights'' it is possible to fit a PARAFAC')
		disp(' model in a weighted least squares sense')
		disp(' ')
		disp(' The input is an array of the same size as the input data holding individual')
		disp(' weights for changing the loss function from a least squares to weighted least')
		disp(' squares one. Instead of minimizing the frobenius norm ||x-M|| where M is the PARAFAC')
		disp(' model, the norm ||(x-M).*weights|| is minimized. The algorithm used for')
		disp(' weighted regression is a majorization step according to Kiers, Psychometrika,')
		disp(' 1997, 62, 251-266, which has the advantage of being computationally cheap and ')
		disp(' also handles situations where ordinary weighted least squares regression does')
		disp(' not work.')
	end
	return
end



% INITIALIZATION OF ALGORITHM
Show = 100; % How often are fit values shown
xsize = size(x);
order = ndims(x);
if (nargin < 3 | ~strcmp(class(scl),'cell'))
	scl = cell(1,order);
end
for ii=1:order   %added 3/1/00 nbg
	if isempty(scl{ii})|(length(scl{ii})~=xsize(ii))
		scl{ii} = 1:xsize(ii);
	end
end
standardtol = [1e-6 1e-6 10000];
if (nargin < 4 | length(tol)==0)
	tol = standardtol;
else % Changed error rb 11/03/00
	tol(find(tol==0)) = standardtol(find(tol==0));
end
if nargin < 6
	options = [1 1];
end
if length(options)<2
	options = [options ones(1,2-length(options))];
end
ls = options(1); %For line-searh
DumpToScreen = options(2);
if nargin < 7|length(plots)==0
	plots = 1;
end
if nargin < 8
	constraints = zeros(1,order);
end
if length(constraints)~=order
	constraints = zeros(1,order);
end

if nargin < 9 
	DoWeight = 0;
else
	if length(size(weights))~=order 
		DoWeight = 0;
	elseif all(size(weights)==xsize)
		DoWeight = 1;
		WMax     = max(abs(weights(:)));
		W2       = weights.*weights;
	else
		DoWeight = 0;
	end
end

%Iterative preproc
iteroffsets = 0;
if exist('iteroff')==1
	if iscell(iteroff)
		iteroffsets = 1;
	elseif isnumeric(iteroff) % Then convert to cell
		hh = iteroff;
		clear iteroff
		iteroff{1} = hh;
		iteroffsets = 1;
	elseif ischar(iteroff)
		if strcmp(lower(iteroff(1:2)),'he')
			clear iteroff
			button = questdlg('Do you want help to incorporate offsets?','Offsets','Yes','No','Cancel','Yes');
			if strcmp(button,'Cancel')
				break
			elseif strcmp(button,'No')
				iteroffsets = 0;
			else
				Continue = 1;
				while Continue
					if Continue==1
						questdlg('For each type of offset, you input the number of the modes across which the offsets are constant. E.g. for a three-way array where you want ordinary offsets that are constant across the first mode (~average of each column); then you input the offset as [1], meaning the the offset is constant across the first mode. In this case the offsets will therefore be a matrix of dimension J x K (J second mode dimension, K third mode dimension). If you want an offset that is constant across all three modes, you input [1 2 3], implying that one offset is estimated','Offsets - how to do it','OK','OK');
					end
					prompt={['Enter the definition of offset # ',num2str(Continue)]};
					def={'[1]'};
					dlgTitle='Offset definition';
					lineNo=1;
					% Ask for offset
					Def=inputdlg(prompt,dlgTitle,lineNo,def);
					% If cancel is hit instead of OK
					if isempty(Def)
						hh=questdlg('Offsets across the first mode will be assumed','Drop offsets?','OK','Stop PARAFAC','OK');
						if strcmp(hh,'Stop PARAFAC')
							error(' PARAFAC stopped')
						else
							Def{1} = '[1]';
						end
					end
					% Verify the offset and show how it's implemented
					hh = eval(Def{1});
					txt = ['[',num2str(hh(1))];
					for i = 2:length(hh);
						txt = [txt,' ',num2str(hh(i))];
					end
					txt = [txt,']'];
					prompt={['Offset # ',num2str(Continue),' is defined as: iteroff{',num2str(Continue),'} = ',txt,';']};
					dlgTitle='Learn how to use the offset input <iteroff>';
					% Ask for offset
					questdlg(prompt,dlgTitle,'OK','Cancel','OK');
					iteroff{Continue} = eval(Def{1});
					Continue = Continue + 1;
					answer=questdlg('Do you want to incorporate other offsets?','More?','Yes','No','Yes');
					if strcmp(answer,'No')
						Continue = 0;
					end
				end
				iteroffsets = 1;
			end
		end
	end
end
if iteroffsets
	XTrilin = zeros(size(x));
	%   options(1) = 0; % No line-search can be performed as iterative offsets are not currently supported in the line-search
end

FeasibilityProblems = 0; % for indicating if the algorithm hasn't yet reached a feasible solution due some constraints
%Define in which mode the scale should go - default is last unless it's fixed')
if all(constraints==-1)
	error(' All modes are fixed, hence no fitting problem')
else
	if constraints(order)==-1 % Since last mode is fixed, scale of the components cannot go there. Instead use the first non-fixed mode
		ScaleMode = min(find(constraints~=-1));
	else
		ScaleMode = order;
	end
end

% CHECK FOR MISSING
if any(isinf(x(:)))|any(isnan(x(:)))
	MissId    = sparse(find(isinf(x)|isnan(x)));
	Missing   = 1;
else
	Missing = 0;
end


% INITIALIZE LOADINGS (if not already given)
if (nargin < 5 | ~(strcmp(class(x0),'cell') | strcmp(class(x0),'struct')))  % Old loadings not given
	x0 = cell(1,order);
	if all(xsize>nocomp) & ~Missing% Use atld
		x0=atld(x,nocomp,0);
		InitString = ' Using fast approximation for initialization';
	elseif order == 3&~Missing&~all(xsize<nocomp)
		% Initialize with TLD estimates
		m = tld(x,nocomp,0,0);
		x0 = m.loads;
		InitString = ' Using direct trilinear decomposition for initialization';
	else
		for j = 1:order
			x0{1,j} = rand(xsize(j),nocomp);
		end
		InitString = ' Using random values for initialization';
	end
else % Old loadings given
	if strcmp(class(x0),'cell')
		x0 = x0;
		InitString = ' Using old values for initialization';
	elseif strcmp(class(x0),'struct')
		if ~strcmp(x0.name,'PARAFAC')
			error(' Input x0 is model that is not a PARAFAC model (name in structure should be PARAFAC)')
		else
			% Use prior model given in x0 to extract loadings and possibly offsets
			if isfield(x0,'Offsets')  % If offsets are given in prior model, then include these
				iteroffsets = 1;
				iteroff = x0.Offsets.Def;
				for i=1:length(iteroff);
					Offset{i}.OffsetsReduced = x0.Offsets.Parameters{1};
					Offset{i}.Offsets = x0.Offsets.Model{1};
				end
				x0 = x0.loads;
				XTrilin = outerm(x0);%zeros(size(x));
				InitString = ' Using old values for initialization (including offsets)';
			else
				x0 = x0.loads;
				InitString = ' Using old values for initialization';
			end
		end
	elseif ~strcmp(class(x0),'cell')
		error('Initial estimate x0 not a cell array')
	end
end


% CHECK FOR INFEASIBLE INITIAL SOLUTIONS
if any(constraints)==1|any(constraints==2) % Nonnegativity required => make sure initial values are positive
	for j=1:order
		if constraints(j)==1|constraints(j)==2
			if any(x0{j}(:)<0) % Negative values occur
				x0{j} = sign(sum(sum(x0{j})))*x0{j}; %'Turn' it around so mainly positive
				i=find(x0{j}<0);
				x0{j}(i)=-x0{j}(i);
			end
		end
	end
end
if any(constraints==3) % Orth required
	for j=1:order
		if constraints(j)==3
			if size(x,j)<nocomp
				error([' Orthogonality cannot be applied in a mode (mode ',num2str(j),') where the dimension is less than the number of factors'])
			else
				x0{j} = orth(x0{j});
			end
		end
	end
end

if Missing % exchange missing elements with model estimates
	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);
	XTrilin = xest; % Multilinear part to subtract from data before calculating offsets
	if iteroffsets & exist('Offset') % it doesn't exist unless prior model is input
		for i=1:length(iteroff) % estimate each set of offsets one at a time
			xest = xest + Offset{i}.Offsets;
		end
	end
	x(MissId)=xest(MissId);
end
% Calculate total sum of squares in data
if DoWeight
	xsq = (x.*weights).^2;
else
	xsq = x.^2;
end
if Missing
	xsq(MissId)=0;  
end
tssq = sum(xsq(:));

% Initialize the unfolded matrices
xuf = cell(1,order);
xuflo = cell(1,order);
xufsize = zeros(1,order);
for i = 1:order
	xuf{i} = unfoldmw(x,i)';
	xufsize(i) = prod(xsize)/xsize(i);
	% Old init of loads;xuflo{i} = zeros(prod(xsize)/xsize(i),nocomp);
	for j = 1:nocomp
		% Old version if i == 1,         mwvect = x0{2}(:,j);,         for k = 3:order              mwvect = mwvect*x0{k}(:,j)';            mwvect = mwvect(:);         end      else         mwvect = x0{1}(:,j);         for k = 2:order            if k ~= i               mwvect = mwvect*x0{k}(:,j)';               mwvect = mwvect(:);            end         end      end      xuflo{i}(:,j) = mwvect;   end
		xuflo{j}=outerm(x0,j,1);
	end
end

% Initialize other variables needed in the ALS
oldx0 = x0;
searchdir = x0;
iter = 0;
flag = 0;
oldests = zeros(prod(xsize)*nocomp,1);
abschange = 0;
relchange = 0;


% Show algorithmic settings etc.
if DumpToScreen
	disp(' ')
	disp(' Fitting PARAFAC ...') 
	txt=[];
	for i=1:order-1
		txt=[txt num2str(xsize(i)) ' x '];
	end
	txt=[txt num2str(xsize(order))];
	disp([' Input: ',num2str(order),'-way ',txt, ' array'])
	
	disp([' A ',num2str(nocomp),'-component model will be fitted'])
	for i=1:order
		if constraints(i)==0
			disp([' No constraints on mode ',num2str(i)])
		elseif constraints(i)==1
			disp([' Nonnegativity on mode ',num2str(i)])
		elseif constraints(i)==2
			disp([' Unimodality on mode ',num2str(i)])
		elseif constraints(i)==3
			disp([' Orthogonality on mode ',num2str(i)])
		elseif constraints(i)==-1
			disp([' Fixed loadings in mode ',num2str(i)])
		end
	end
	
	disp([' Convergence criteria:']) 
	disp([' Relative change in fit : ',num2str(tol(1))]) 
	disp([' Absolute change in fit : ',num2str(tol(2))]) 
	disp([' Maximum iterations     : ',num2str(tol(3))]) 
	
	if Missing
		disp([' ', num2str(100*(length(MissId)/prod(xsize))),'% missing values']);
	else
		disp(' No missing values')
	end
	
	if DoWeight
		disp(' Weighted optimization will be performed using input weights')
	end
	
	%Iterative preproc
	if iteroffsets
		disp(' Iteratively calculated offsets will be used')
		disp(' Line-search has been turned off')
	else
		%disp(' No iterative offsets')
	end
	
	disp(InitString)
end

% Start the ALS
while flag == 0;
	iter = iter+1;
	% Loop over each of the order to estimate
	
	if DoWeight & iter > 1% If Weighted regression is to be used, do majorization to make a transformed data array to be fitted in a least squares sense
		out = reshape(xest,xsize) + (WMax^(-2)*W2).*(x - reshape(xest,xsize));
		for i = 1:order
			xuf{i} = unfoldmw(out,i)';
		end
		clear out
	end
	
	%Iterative preproc
	if iteroffsets
		if iter>1
			OldOffset = Offset;
		end
		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) - Xtrilin;
		else
			xprep = x - XTrilin;
		end