gemanova.m

进行方差分析的MATLAB源码

% Together B1 and BI can be concatenated to give the solution to
% problem ULSR for any modloc position AS long as we do not pay
% attention to the element of x at this position


All=zeros(I,I+2);
All(1:I,3:I+2)=B1;
All(1:I,1:I)=All(1:I,1:I)+BI;
All=All(:,2:I+1);
Allmin=All;
Allmax=All;
% All(:,i) holds the ULSR solution for modloc = i, disregarding x(i),


iii=find(x>=max(All)');
b=All(:,iii(1));
b(iii(1))=x(iii(1));
Bestfit=sum((b-x).^2);
MaxML=iii(1);
for ii=2:length(iii)
	this=All(:,iii(ii));
	this(iii(ii))=x(iii(ii));
	thisfit=sum((this-x).^2);
	if thisfit<Bestfit
		b=this;
		Bestfit=thisfit;
		MaxML=iii(ii);
	end
end

if xmin<0
	b=b+xmin;
end


% Impose nonnegativity
if NonNeg==1
	if any(b<0)
		id=find(b<0);
		% Note that changing the negative values to zero does not affect the
		% solution with respect to nonnegative parameters and position of the
		% maximum.
		b(id)=zeros(size(id))+0;
	end
end



function [b,B,AllBs]=monreg(x);

% Monotonic regression according
% to J. B. Kruskal 64
%
% b     = min|x-b| subject to monotonic increase
% B     = b, but condensed
% AllBs = All monotonic regressions, i.e. AllBs(1:i,i) is the 
%         monotonic regression of x(1:i)
%
%
% Copyright 1997
%
% Rasmus Bro
% Royal Veterinary & Agricultural University
% Denmark
% [e-mail]
%


I=length(x);
if size(x,2)==2
	B=x;
else
	B=[x(:) ones(I,1)];
end

AllBs=zeros(I,I);
AllBs(1,1)=x(1);
i=1;
while i<size(B,1)
	if B(i,1)>B(min(I,i+1),1)
		summ=B(i,2)+B(i+1,2);
		B=[B(1:i-1,:);[(B(i,1)*B(i,2)+B(i+1,1)*B(i+1,2))/(summ) summ];B(i+2:size(B,1),:)];
		OK=1;
		while OK
			if B(i,1)<B(max(1,i-1),1)
				summ=B(i,2)+B(i-1,2);
				B=[B(1:i-2,:);[(B(i,1)*B(i,2)+B(i-1,1)*B(i-1,2))/(summ) summ];B(i+1:size(B,1),:)];
				i=max(1,i-1);
			else
				OK=0;
			end
		end
		bInterim=[];
		for i2=1:i
			bInterim=[bInterim;zeros(B(i2,2),1)+B(i2,1)];
		end
		No=sum(B(1:i,2));
		AllBs(1:No,No)=bInterim;
	else
		i=i+1;
		bInterim=[];
		for i2=1:i
			bInterim=[bInterim;zeros(B(i2,2),1)+B(i2,1)];
		end
		No=sum(B(1:i,2));
		AllBs(1:No,No)=bInterim;
	end
end

b=[];
for i=1:size(B,1)
	b=[b;zeros(B(i,2),1)+B(i,1)];
end


function [x,w] = CrossProdFastNnls(XtX,Xty,tol)
%NNLS	Non-negative least-squares.
%	b = CrossProdFastNnls(XtX,Xty) returns the vector b that solves X*b = y
%	in a least squares sense, subject to b >= 0, given the inputs
%       XtX = X'*X and Xty = X'*y.
%
%	[b,w] = fastnnls(XtX,Xty) also returns dual vector w where
%	w(i) < 0 where b(i) = 0 and w(i) = 0 where b(i) > 0.
%	L. Shure 5-8-87 Copyright (c) 1984-94 by The MathWorks, Inc.
%
%  Revised by:
%	Copyright
%	Rasmus Bro 1995
%	Denmark
%	E-mail [e-mail]
%  According to Bro & de Jong, J. Chemom, 1997, 11, 393-401

% initialize variables


if nargin < 3
	tol = 10*eps*norm(XtX,1)*max(size(XtX));
end
[m,n] = size(XtX);
P = zeros(1,n);
Z = 1:n;
x = P';
ZZ=Z;
w = Xty-XtX*x;

% set up iteration criterion
iter = 0;
itmax = 30*n;

% outer loop to put variables into set to hold positive coefficients
while any(Z) & any(w(ZZ) > tol)
	[wt,t] = max(w(ZZ));
	t = ZZ(t);
	P(1,t) = t;
	Z(t) = 0;
	PP = find(P);
	ZZ = find(Z);
	nzz = size(ZZ);
	z(PP')=(Xty(PP)'/XtX(PP,PP)');
	z(ZZ) = zeros(nzz(2),nzz(1))';
	z=z(:);
	% inner loop to remove elements from the positive set which no longer belong
	
	while any((z(PP) <= tol)) & iter < itmax
		
		iter = iter + 1;
		QQ = find((z <= tol) & P');
		alpha = min(x(QQ)./(x(QQ) - z(QQ)));
		x = x + alpha*(z - x);
		ij = find(abs(x) < tol & P' ~= 0);
		Z(ij)=ij';
		P(ij)=zeros(1,max(size(ij)));
		PP = find(P);
		ZZ = find(Z);
		nzz = size(ZZ);
		z(PP)=(Xty(PP)'/XtX(PP,PP)');
		z(ZZ) = zeros(nzz(2),nzz(1));
		z=z(:);
	end
	x = z;
	w = Xty-XtX*x;
end

x=x(:);





function loads=atld(X,F,Show);



if nargin<3
	Show = 0;
end

xsize = size(X);
order = ndims(X);
RealData = all(isreal(X(:)));

% Initialize with random numbers
for i = 1:order;
	if xsize(i)>=F
		loads{i} = orth(rand(xsize(i),F));
	else
		loads{i} = rand(xsize(i),F);
	end
	invloads{i} = pinv(loads{i});
end

model = outerm(loads);
fit=sum(abs((X(:)-model(:)).^2));
oldfit=2*fit;
maxit=30;
crit=1e-6;
it=0;

while abs(fit-oldfit)/oldfit>crit&it<maxit
	it=it+1;
	oldfit=fit;
   
   
	% Normalize loadings   
	for mode = 2:order
		scale = sqrt(abs(sum(loads{mode}.^2)));
		loads{1}    = loads{1}*diag(scale);
		loads{mode} = loads{mode}*diag(scale.^-1);
		% Scale occasionally to mainly positiviy
		if RealData & (it==1|it==2|rem(it,1)==0)
			scale = sign(sum(loads{mode}));
			loads{1}    = loads{1}*diag(scale);
         loads{mode} = loads{mode}*diag(scale);
         if any(isnan(loads{mode}))
            loads{mode} = rand(size(loads{mode}));
         end
			invloads{mode} = pinv(loads{mode});
		end
	end
	
	if it == 1
		delta = linspace(1,F^(order-1),F);
	end
	
	% Compute new loadings for all modes
	for mode = 1:order;
		
		xprod = X;
		% Multiply pseudo-inverse loadings in all but the mode being estimated
		if mode == 1
			for mulmode = 2:order
				xprod = ntimes(xprod,invloads{mulmode},2,2);
			end
		else
			for j = 1:mode-1
				xprod = ntimes(xprod,invloads{j},1,2);
			end
			for j = mode+1:order
				xprod = ntimes(xprod,invloads{j},2,2);
			end
		end
		
		% Extract first mode loadings from product
		loads{mode} = xprod(:,delta);
		invloads{mode} = pinv(loads{mode});
	end
	
	
	model = outerm(loads);
	fit=sum((X(:)-model(:)).^2);
end
% Normalize loadings   
for mode = 2:order
	scale = sqrt(sum(loads{mode}.^2));
	loads{1}    = loads{1}*diag(scale);
	loads{mode} = loads{mode}*diag(scale.^-1);
end

if Show
	disp(' ')
	disp('    Iteration    sum-sq residuals')
	disp(' ')
	fprintf(' %9.0f       %12.10f    \n',it,fit);
end




function product = ntimes(X,Y,modeX,modeY);


%NTIMES Array multiplication
% X*Y is the array/matrix product of X and Y. These are multiplied across the 
% modeX mode/dimension of X and modeY mode/dimension of Y. The number of levels of X in modeX
% must equal the number of levels in modeY of Y.
%
% The product will be an array of order two less than the sum of the orders of A and B
% and thus works as a straightforward extension of the matrix product. The order
% of the modes in the product are such that the X-modes come first and then the
% Y modes. 
%
% E.g. if X is IxJ and Y is KxL and I equals L, then 
% ntimes(X,Y,1,2) will yield a JxK matrix equivalent to the matrix product
% X'*Y'. If X is IxJxK and Y is LxMxNxO with J = N then 
% ntimes(X,Y,2,3) yields a product of size IxKxLxMxO
% 
%I/O: product = NTIMES(X,Y,modeX,modeY) 
% 
%See also TIMES,MTIMES.

%Copyright Eigenvector Research Inc./Rasmus Bro, 2000
%Rasmus Bro, August 20, 2000

orderX = ndims(X);
orderY = ndims(Y);
xsize  = size(X);
ysize  = size(Y);

X = permute(X,[modeX 1:modeX-1 modeX+1:orderX]);
Y = permute(Y,[modeY 1:modeY-1 modeY+1:orderY]);
xsize2  = size(X);
ysize2  = size(Y);

if size(X,1)~=size(Y,1)
	error(' The number of levels must be the same in the mode across which multiplications is performed')
end

%multiply the matricized arrays
product = reshape(X,xsize2(1),prod(xsize2(2:end)))'*reshape(Y,ysize2(1),prod(ysize2(2:end)));
%reshape to three-way structure
product = reshape(product,[xsize2(2:end) ysize2(2:end)]);



function mwauf = unfoldmw(mwa,order,meth)
%UNFOLDMW Unfolds multiway arrays along specified order
% The inputs are the multiway array to be unfolded (mwa),
% and the dimension number along which to perform the
% unfolding (order). The output is the unfolded array (mwauf).
% This function is used in the development of PARAFAC models
% in the alternating least squares steps.
%
%I/O: mwauf = unfoldmw(mwa,order);
%
%See also: MPCA, OUTER, OUTERM, PARAFAC, TLD

%Copyright Eigenvector Research, Inc. 1998
%bmw
%rb August 2000, speeded up by factor 20


mwasize = size(mwa);
ord = length(mwasize);
mwauf=permute(mwa,[order 1:order-1 order+1:ord]);
mwauf=reshape(mwauf,mwasize(order),prod(mwasize([1:order-1 order+1:ord])));

%Old version
%mwasize = size(mwa);
%ms = mwasize(order);
%po = prod(mwasize);
%ns = po/ms;
%if order ~= 1
%   pod = prod(mwasize(1:order-1));
%end
%mwauf = zeros(ms,ns);
%for i = 1:ms
%   if order == 1
%      mwauf(i,:) = mwa(i:ms:po);
%   else
%      inds = zeros(1,ns); k = 1; fi = (i-1)*pod + 1;
%      for j = 1:ns/pod
%         inds(k:k+pod-1) = fi:fi+pod-1;
%         fi = fi + ms*pod;
%         k = k + pod;
%      end
%      mwauf(i,:) = mwa(inds);
%   end
%end%
%
%end

function mwa = outerm(facts,lo,vect)
%OUTERM Outer product of any number of vectors with multiple factors
%  The input to outerm is a 1 by n cell array (facts), where each cell
%  contains the factors for one of the ways, or orders, with each
%  of the factors being a column in the matrix. Optional inputs
%  are the number of an order to leave out (lo) in the formation
%  of the product, and a flag (vect) which causes the function
%  to not sum and reshape the final factors when set to 1. (This option
%  is used in the alternating least squares steps in PARAFAC.) 
%  The output is the multiway array resulting from multiplying the
%  factors together(mwa), or the strung out individual factors.
%
%I/O: mwa = outerm(facts,lo,vect);
%
%See also: OUTER, PARAFAC, TLD

%Copyright Eigenvector Research, Inc. 1998
%bmw

if nargin < 2
  lo = 0;
end
if nargin < 3
  vect = 0;
end
order = length(facts);
if lo == 0
  mwasize = zeros(1,order);
else
  mwasize = zeros(1,order-1);
end
k = 0;
for i = 1:order
  [m,n] = size(facts{i});
  if i ~= lo
    k = k + 1;
    mwasize(k) = m;
  end
  if i > 1
    if nofac ~= n
	  error('All orders must have the same number of factors')
	end
  else
    nofac = n;
  end
end
mwa = zeros(prod(mwasize),nofac);

for j = 1:nofac
  if lo ~= 1
    mwvect = facts{1}(:,j);
    for i = 2:order
	  if lo ~= i
        %mwvect = kron(facts{i}(:,j),mwvect);
		mwvect = mwvect*facts{i}(:,j)';
		mwvect = mwvect(:);
	  end
    end
  elseif lo == 1
    mwvect = facts{2}(:,j);
	for i = 3:order
      %mwvect = kron(facts{i}(:,j),mwvect);
	  mwvect = mwvect*facts{i}(:,j)';
	  mwvect = mwvect(:);
	end
  end
  mwa(:,j) = mwvect;
end
% If vect isn't one, sum up the results of the factors and reshape
if vect ~= 1
  mwa = sum(mwa,2);
  mwa = reshape(mwa,mwasize);
end