gemanova.m

进行方差分析的MATLAB源码

function model = gemanova(X,F,Fix,scl,cross,show)

%GEMANOVA - GEneralized Multiplicative ANOVA
% Fits a GEMANOVA model to an N-way array of responses 
% X with F effects. 
%
% INPUTS
% The third input, Fix defines the nature 
% of the effects. Fix is a binary matrix of size 
% N x F and the f'th column defines the f'th effect. 
% If the column contain zeros only, the effect is
% a multiplicative effect of all modes.
%
% Ex.: For a three-way array of size I x J x K, 
%
% Fix(:,f) = [0 0 0]' corresponds to a_i*b_j*c_k
% Fix(:,f) = [1 0 0]' corresponds to b_j*c_k
% Fix(:,f) = [1 0 1]' corresponds to b_j
% Fix(:,f) = [1 1 1]' corresponds to m (an overall constant)
% 
% Thus a two-effect model: x_ijk = a_i*b_j*c_k + c_k
% can be specified as
%
% model = gemanova(X,2,[0 0 0;1 1 0]',scl);
%
% scl is an optional input of a cell array of vectors for 
% plotting effects against (set scl=[] for none)
%
% A fifth optional input, cross, performs cross-validation
% if set to one. 
% 
% OUPUTS
% The output is a structured array model which holds
%     xname: name of the original workspace input variable
%      name: type of model, always 'GEMANOVA'
%      date: model creation date stamp
%      time: model creation time stamp
%      size: size of the original input array
%  neffects: number of effects estimated
%       ssq: residual sum of squares
%   effects: 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 
%
% If cross-validation is performed, the additional output
% is a structure Cross which holds the reproduced 
% array obtained by leave-one-element-out cross-validation 
% (as the elements are assumed to be independent in ANOVA) 
% as well as the RMSECV and the correlation between the 
% responses and the reproduced array. Note that MANOVA can be
% performed directly with gemanova but the cross-validation
% assumes independent elements.
%
% The algorithm handles missing elements if these are 
% set to NaN of Inf. 
%
% Please refer to this m-file through
%
%      Rasmus Bro & Marianne Jakobsen, Exploring complex 
%      interactions in designed data using GEMANOVA. Color 
%      changes in fresh beef during storage, Journal of 
%      Chemometrics, 2001, Submitted.
%
% I/O: model = gemanova(X,2,[0 0 0;1 1 0]',scl,Cross);

%rb February, 2001

id = find(isinf(X));
X(id) = NaN; % Because the following algorithms only handle NaN

if (nargin < 4 | ~strcmp(class(scl),'cell'))
    scl = cell(1,length(size(X)));
end
if nargin < 5
    cross = 0;
end
if nargin < 6
    show = 1;
end
if any(sum(Fix)==length(size(X)));
    % Means that there is a constant across all modes
    % Then an extra mode of dimension one has to be added
    % in the end, the mode will be removed again
    x = zeros([1 size(X)]);
    x(1,:) = X(:)';
    X = x;
    clear x;
    Fix = [ones(1,F);Fix];
    i = find(sum(Fix)==length(size(X)));
    Fix(1,i)=0;
    ExtraMode = 1;
else
    ExtraMode = 0;
end

    



if cross==1
    % Find overall model
    disp(' ')
    disp(' CROSS-VALIDATING THE SOUGHT MODEL')
    disp('  Fitting the global model ...')
    model = gemanova(X,F,Fix,scl);
    disp(' ')
    disp('  Cross-validating ...')
    Cross = gemcross(X,F,Fix,scl,model.nparam);
    model.Cross = Cross;
    
    
else % Fit model
    
    if show
        disp(' ')
        disp(' Gradually adjusting to correct model')
    end
    model = parafacweight(X,F,Fix,show);
    
    if show
        disp(' Refining to exact solution')
    end
    [m,fit]=parafacones(X,F,Fix,5000,model,show);
    
    % fix added mode if necessary
    if ExtraMode
        Fix = Fix(2:end,:);
        clear effect
        for i = 1:length(m)-1
            effect{i} = m{i+1};
        end
        effect{end}=effect{end}*diag(m{1});
        m = effect;
        ExtraMode = 0;
        X = squeeze(X);
        if show
            disp(' Offsets constant across all modes are incorporated into the last mode parameters')
        end
        
    end
    
    
    % OUTPUT MODEL
    model = struct('xname',inputname(1),'name','GEMANOVA','date',date,'time',clock,...
        'size',size(X),'neffects',F);
    model.effects = m;
    model.ssq = fit;
    model.scl = scl;
    model.res = X-outerm(model.effects);
    nparam = 0;
    for i = 1:size(Fix,1)
        nparam = nparam + sum(Fix(i,:)==0)*size(X,i);
    end
    model.nparam = nparam;
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% AUXILIARY FILES%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Cross = gemcross(X,F,Fix,scl,nparam);

% Do cross-validation
Idx = find(~isnan(X));
Pred = NaN*ones(size(X));
disp([' ',num2str(length(Idx)),' segments']);
fprintf([' '])
FIT = [];
for i = 1:length(Idx)
    x = X;
    x(Idx(i)) = NaN;
    crossmodel = gemanova(x,F,Fix,scl,0,0);
    crossmod = outerm(crossmodel.effects);
    Pred(Idx(i))=crossmod(Idx(i));
    e = X(find(~isnan(x)))-crossmod(find(~isnan(x)));
    fit = sum(e(:).^2);
    FIT = [FIT;fit]; % Save this to check that there are no strange fit values.
    if i~=length(Idx)
        fprintf([num2str(i),','])
    else
        fprintf(num2str(i))
    end
%    disp([' Crossval segment ',num2str(i),' of ',num2str(length(Idx)),' - ssq: ',num2str(sum(e(:).^2))]);
end
Cross.Predict = Pred;
Cross.RMSECV = std(X(Idx)-Pred(Idx));
e = (X(Idx)-Pred(Idx)).^2;
Cross.RMSECVcorrect = sqrt( sum(e)/(length(Idx)-nparam));
r = corrcoef([vec(X(Idx)) vec(Pred(Idx))]);
Cross.corr = r(2,1);
Cross.IndividualSsq = FIT;




function model = parafacweight(X,F,Fix,show)

% Parafac imposing fixed elements with ones gradually through weighting
% The end result will not have exact ones but will be a good initial 
% starting point for an exact algoritm


crit = 1e-10;
showfit = pi; % Don't show anything
showfitinit = pi;
numsplitgradient = 6;
maxit = 500;
if nargin<4
    show = 1;
end
order = length(size(X));
it = 0;
if size(Fix,1)~=order & size(Fix,2)~=F
    error(' Fix is not correctly given')
end
DimX = size(X);


% Initialize with parafac
fit = sum(X(find(~isnan(X))).^2);
fitold = fit*2;
m = parafac(X,F,[],[],[],[1 0],0);
for rep = 1:3 % See if randomly started is better
    x0 = m.loads;
    for i = 1:length(x0)
        x0{i} = rand(size(x0{i}));
    end
    m2 = parafac(X,F,[],[],x0,[1 0],0);
    if m2.ssq(2)<m.ssq(2)
        m = m2;
    end
end
        
loads = m.loads;

% Find gradually increasing weights
ssqX = fit;
Weights = linspace(0,ssqX/300,numsplitgradient);

for w = 1:numsplitgradient
    fitold = fit*10;
    it = 0;
    while abs(fit-fitold)/fitold>crit & it < maxit
        
        it = it+1;
        fitold = fit;
        
        % Adjust scales on different loads to avoid numerical problems
        scales = ones(F,1);
        for fac = 1:F
            for ord = 1:order
                if ~Fix(ord,fac)
                    scales(fac) = scales(fac)*norm(loads{ord}(:,fac));
                end
            end
        end
        for fac = 1:F
            NumbNonFix = length(find(~Fix(:,fac)));
            scales(fac) = scales(fac)^(1/NumbNonFix);
            for ord = 1:order
                if ~Fix(ord,fac)
                    loads{ord}(:,fac) = scales(fac)*loads{ord}(:,fac)/norm(loads{ord}(:,fac));
                end
            end
        end
        
        for fac = 1:F
            for ord = 1:order
                % find update for ord'th mode of fac'th factor
                z = outerm(loads,ord,1);
                permX = permute(X,[ord 1:ord-1 ord+1:order]);
                permXunf = reshape(permX,DimX(ord),prod(DimX([1:ord-1 ord+1:end])));
                % handle missing
                for p = 1:DimX(ord)
                    id = find(~isnan(permXunf(p,:)));
                    xx= permXunf(p,id)';
                    zz = z(id,:);
                    if Weights(w)&any(Fix(ord,:))
                        xx = [xx;ones(length(find(Fix(ord,:))),1)*Weights(w)];
                        ffix = eye(F);
                        ffix(find(~Fix(ord,:)),:)=[];
                        zz  = [zz;ffix*Weights(w)];
                    end
                    ll = xx'*pinv(zz)';
                    if sum(abs(ll))<100*eps; % almost zero
                        loads{ord}(p,:) = .9*loads{ord}(p,:)+.1*ll;
                    else
                        loads{ord}(p,:) = ll;
                    end
                end
            end
            
        end
        
        
        fit = pffit(X,loads);
        if show & rem(it,showfit)==0 & it > 1
            disp([' Fit of GEMANOVA model = ',num2str(fit),' after ',num2str(it),' it. in final round & Weight = ',num2str(Weights(w)),' (',num2str(w),'/',num2str(numsplitgradient),') '])
        end
    end
end

model = loads;
if show
    disp([' Final fit of GEMANOVA model = ',num2str(fit),' after ',num2str(it),' it. & Weight = ',num2str(Weights(w)),' (',num2str(w),'/',num2str(numsplitgradient),')'])
end
%%%%%%%%%%%%%%%%%%% END OF PROGRAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function ssq = pffit(X,loads);
M = outerm(loads);
E = X-M;
ssq = sum(E(find(~isnan(E))).^2);


function product = khatri(loads,factor);
product = [];
for fac = 1:size(loads{1},2)
    if fac~=factor
        Z = kron(loads{end}(:,fac),loads{end-1}(:,fac));
        for j = length(loads)-2:-1:1
            Z = kron(Z,loads{j}(:,fac));
        end
        product = [product Z];
    end
end


function [m,fit]=parafacones(X,F,Fix,maxit,loads,show)

% Fit exact ls parafac with exact ones

crit = 1e-11;
showfit = pi;
showfitinit = pi;
if nargin < 4
    maxit = 1000;
end
if nargin<6
    show = 1;
end
order = length(size(X));
it = 0;
if size(Fix,1)~=order & size(Fix,2)~=F
    error(' Fix is not correctly given')
end
DimX = size(X);


fit = sum(X(find(~isnan(X))).^2);
fitold = fit*2;
while abs(fit-fitold)/fitold>crit & it < maxit
    
    it = it+1;
    fitold = fit;
    
    if it >5 & rem(it,10)==0 %& 5==3% Do line-search
        options = optimset;
        options = optimset(options,'Display','off');
        [alpha,fithere,exitf] = fminbnd('pffitalpha',1,100,options,X,loads,LoadingsOld);
        if fithere<fit
            for ord = 1:order
                loads{ord} = loads{ord}+alpha*(loads{ord}-LoadingsOld{ord});
            end
            fit = fithere;
        end
        if show & rem(it,showfit)==0
            disp([' Fit after line-search = ',num2str(fit),' after ',num2str(it),' iterations'])
        end
    end
    
    LoadingsOld = loads;
    for fac = 1:F
        
        % Calculate residual of X for model with all but fac'th factor
        if F>1
            product = khatri(loads,fac);
        else
            product = zeros(size(X(:)));
        end
        if F>2
            sumprod = sum(product');
        else
            sumprod = product;
        end
        resX = X - reshape(sumprod,size(X));
        
        % calculate old loads for fac
        for ord = 1:order;
            oldloads{ord} = loads{ord}(:,fac);
        end
        
        
        for ord = 1:order
            % find update for ord'th mode of fac'th factor
            z = outerm(oldloads,ord,1);
            if ~Fix(ord,fac)
                permresX = permute(resX,[ord 1:ord-1 ord+1:order]);
                permresXunf = reshape(permresX,DimX(ord),prod(DimX([1:ord-1 ord+1:end])));
                % handle missing
                for p = 1:DimX(ord)
                    id = find(~isnan(permresXunf(p,:)));
                    oldloads{ord}(p) = permresXunf(p,id)*pinv(z(id))';
                end
            else
                oldloads{ord}=ones(DimX(ord),1);
            end
        end
        
        % update loads
        for ord = 1:order
            loads{ord}(:,fac) = oldloads{ord};
         end
         
    end
    
    
    fit = pffit(X,loads);
    if show & rem(it,showfit)==0 & it > 1
        disp([' Fit of GEMANOVA model = ',num2str(fit),' after ',num2str(it),' iterations'])
     end
end

% Scale all but last (non-fixed) mode to unit length
for f = 1:F
    for o = 1:order-1
        if ~Fix(o,f) % Not fixed to ones, hence should be normalized
           if any(~Fix(o+1:end,f)) % Only normalize if there is a 'later' unfixed mode
              oo = o+1:order;
                fix = find(~Fix(oo,f));
                fix = oo(fix(end));
                loads{fix}(:,f) = loads{fix}(:,f)*norm(loads{o}(:,f));
                loads{o}(:,f) = loads{o}(:,f)/norm(loads{o}(:,f));
            end
        end
    end
end

m = loads;
if show
    disp([' Final fit of GEMANOVA model = ',num2str(fit),' after ',num2str(it),' iterations'])
 end
 
 
 function vX = vec(X);
 
 vX = X(:);