PCA.rar_Principal-component_principal component

function result = mf(X,model,mode,options); %MF Matrix factorization for simplifying PCA and Tucker models. % Takes a PCA or Tucker model and rotates one mode to optimal simplicity. % The input is the data used to fit the original model, the model % structure and the mode to be rotated. The remainder of the model is % counterrotated so that the fit of the model is the same. In PCA % the scores or loadings are counterrotated whereas in Tucker, the core % is counterrotated. % The rotation is done by orthomax where default is to set the parameter % gamma to one (corresponds to varimax) and where setting it to % zero would correspond to quartimax. % Optional structure input (options) can contain any of the fields: % % display [ {'off'} | 'on'] Governs screen display. % gamma between 0 and 1. One is default. % %I/O: result = mf(X,model,mode,options); % %See also: PCA, TUCKER % Copyright � Eigenvector Research, Inc. 1991-2006 % Licensee shall not re-compile, translate or convert "M-files" contained % in PLS_Toolbox for use with any software other than MATLAB�, without % written permission from Eigenvector Research, Inc. if nargin == 0; X = 'io'; end if ischar(X); options = []; options.plots = 'on'; options.display = 'on'; options.preprocessing = {[]}; options.confidencelimit = 0.95; options.gamma = 1; if nargout==0; evriio(mfilename,X,options); else; result = evriio(mfilename,X,options); end return; end if nargin<3 error(' You must define the mode in which you want to rotate the components in MF') end if nargin < 4 | isempty(options); options = mf('options'); else options = reconopts(options,'mf'); end if options.gamma>1 warning('Gamma can at most be one') options.gamma = 1; elseif options.gamma<0 warning('Gamma must be at least zero') options.gamma = 0; end if isa(X,'double') %convert X to DataSet X = dataset(X); X.name = inputname(1); X.author = 'MCR'; elseif ~isa(X,'dataset') error(['Input X must be class ''double'' or ''dataset''.']) end if ndims(X.data)>2 error(['Input X must contain a 2-way array. Input has ',int2str(ndims(varargin{1}.data)),' modes.']) end %Handle Preprocessing if isempty(options.preprocessing); options.preprocessing = {[]}; %reinterpet as empty cell end if ~isa(options.preprocessing,'cell'); options.preprocessing = {options.preprocessing}; %insert into cell end if mdcheck(X.data(X.include{1},X.include{2})); if strcmp(options.display,'on'); warning('Missing Data Found - Replacing with "best guess". Results may be affected by this action.'); end [flag,missmap,replaced] = mdcheck(X.data(X.include{1},X.include{2})); X.data(X.include{1},X.include{2}) = replaced; end %preprocessing if ~isempty(options.preprocessing{1}); [X,options.preprocessing{1}] = preprocess('calibrate',options.preprocessing{1},X); end % GO AHEAD if strcmp(model.modeltype,'PCA') [Arot,T,f]=orthomax(model.loads{mode},options.gamma); if mode==1 othermode = 2; else othermode = 1; end B = model.loads{othermode}*pinv(T'); loads{mode,1} = Arot; loads{othermode,1} = B; % % Fit a two-way PARAFAC to obtain a complete structure % o = parafac('options'); % o.constraints = o.constraints(1:2); % for j = 1:2 % o.constraints{j}.fixed.position = ones(size(result.loads{j})); % o.constraints{j}.fixed.value = result.loads{j}; % o.constraints{j}.fixed.weight = -1; % end % o.display = 'off'; % o.waitbar = 'off'; % o.plots = 'off'; % result = parafac(X,size(result.loads{1},2),o); % result.modeltype = 'matrix factorization'; % det = result.detail; % det = rmfield(det,{'stopcrit','critfinal','options','coreconsistency','algo','initialization'}); % result.detail = det; result = modelstruct('mcr'); result.date = date; result.time = clock; result = copydsfields(X,result); result.loads = loads; result.rotationtype = ['ORTHOMAX gamma = ',num2str(options.gamma)]; result.rotationmatrix = T; result.detail.data{1} = X; result.pred{1} = result.loads{1,1}*result.loads{2,1}'; result.detail.res{1} = X.data(:,result.detail.includ{2,1}) - result.pred{1}; result.ssqresiduals{1,1} = result.detail.res{1}.^2; result.ssqresiduals{2,1} = NaN*ones(1,size(X.data,2)); result.ssqresiduals{2,1}(1,result.detail.includ{2,1}) = ... sum(result.ssqresiduals{1,1}(result.detail.includ{1,1},:),1); %var SSQs based on cal samples only result.ssqresiduals{1,1} = sum(result.ssqresiduals{1,1},2); %sample SSQs for ALL samples (excluded vars already gone) %calculate SSQ table %calculate % signal by factor for j=1:size(result.loads{1},2); sig(j) = sum(sum((result.loads{1}(result.detail.includ{1,1},j)*result.loads{2}(:,j)').^2)); end sig = normaliz(sig,[],1); %normalize to 100% uncap = sum(sum(result.detail.res{1}(result.detail.includ{1,1},:).^2))./sum(sum(mncn(X.data(X.include{1},X.include{2})).^2)); if uncap>=1; %happens with really bad fits (e.g. poor constraints or badly processed data) uncap = sum(sum(result.detail.res{1}(result.detail.includ{1,1},:).^2))./sum(sum(X.data(X.include{1},X.include{2}).^2)); end if uncap>=1; uncap = 0; end result.detail.ssq = [[1:size(result.loads{1},2)]' sig'*100 (1-uncap)*sig'*100 cumsum((1-uncap)*sig')*100]; switch options.display case 'on' ssqtable(result) end if length(result.detail.includ{1,1})>size(B,2) result.detail.tsqlim{1,1} = tsqlim(length(result.detail.includ{1,1}),size(B,2),options.confidencelimit*100); else result.detail.tsqlim{1,1} = 0; end try switch lower(options.plots) case 'final' plotloads(result); plotscores(result); case 'on' plotloads(result); plotscores(result); end catch warning(lasterr) end elseif strcmp(lower(model.modeltype),'tucker') error('Not implemented yet') end function [Arot,T,f]=orthomax(A,gamma); %ORHTOMAX % Orthomax rotation of A % based on Kiers (1997), psychometrika, 62, 579-598 % according to Clarkson % Jennrich (1988), psychometrika, 53(2), 251-259 % Remedy of ten Berge (1995), psychometrika, 60, 437-446 not built in % % [Arot,T,f]=orthomax(A,gamma); % % input % A loading-matrix to rotate % gamma orthomax parameter % % output % Arot rotated A % T rotation matrix % Copyright � Eigenvector Research, Inc. 1991-2006 % Licensee shall not re-compile, translate or convert "M-files" contained % in PLS_Toolbox for use with any software other than MATLAB�, without % written permission from Eigenvector Research, Inc. ConvCrit=1e-6; [I,F]=size(A); Arot = A; f= sum(Arot(:).^4) - gamma/I*sum((sum(Arot.^2)).^2); fold=2*f; it=0; k3 = -gamma/I; k4 = 1; while abs(f-fold)/abs(f)>ConvCrit fold=f; it=it+1; for f1=1:F-1 for f2=f1+1:F x=Arot(:,f1); y=Arot(:,f2); a3 = .25*(x'*x -y'*y)^2-(x'*y)^2; a4 = .25*(x.^2'*x.^2 +y.^2'*y.^2)-1.5*(x.^2'*y.^2); b3 = x'*y*(x'*x-y'*y); b4 = x.^3'*y-x'*y.^3; a = k3*a3+k4*a4; b = k3*b3+k4*b4; theta = .25*atan2(b,a); xnew = x*cos(theta)+y*sin(theta); ynew = -x*sin(theta)+y*cos(theta); u = x.*x-y.*y; v = 2*x.*y; u = u-mean(u); v = v-mean(v); a = 2*I*u'*y; b = I*u'*u-I*y'*y; c = (a^2+b^2)^(.5); v = -sign(a)*( (b+c)/2*c )^(.5); Arot(:,f1) = xnew; Arot(:,f2) = ynew; end end f= sum(Arot(:).^4) - gamma/I*sum((sum(Arot.^2)).^2); end; T = pinv(A)*Arot;