pca-test.rar_PCA数据_pca

clc clear load pcadata; figure() plot(part1); title('Data for PCA example, Part 1') xlabel('Time') ylabel('Variable Values') [m,n] = size(part1); [std_data,mean_part1,std_part1] = standardization(part1); figure() plot(std_data); title('Normalized Data for PCA example, Part 1') xlabel('Time') ylabel('Variable Values') [u,s,v] = svd(std_data); score = u*s; loads = v; singular = diag(s); variance = singular.^2; contr = variance/sum(variance); acc_contr = cumsum(variance)/sum(variance); r=sum(singular~=0); temp2 = (1:n)'; ssq = [temp2 singular contr*100 acc_contr*100]; %print variance information disp(' ') disp(' Percent Variance Captured by PCA Model') disp(' ') disp('Principal Eigenvalue % Variance % Variance') disp('Component of Captured Captured') disp(' Number Cov(X) This PC Total') disp('--------- ---------- ---------- ----------') format = ' %3.0f %3.2e %6.2f %6.2f'; mprint = min([20 n m]); for i = 1:mprint tab = sprintf(format,ssq(i,:)); disp(tab) end figure() stem(1:r,contr,'o'),title('Principal Component Contribution') xlabel('Principal Component Sequence Number'),ylabel('Contribution') figure() plot(1:r,acc_contr,'-o'),title('Principal Component Cumulative Contribution') xlabel('Principal Component Sequence Number'),ylabel('Cumulative Contribution') %选择主元个数 r = input('How many principal components do you want to keep? '); if r > n error('Number of PCs must be <= number of variables') elseif r > m error('Number of PCs must be <= number of samples') elseif r < 1 error('Number of PCs must be > 0') elseif isempty(r) error('Number of PCs must be > 0') end score_new = score(:,1:r); load_new = loads(:,1:r); for i = 1:m Q_PCA(i) = std_data(i,:)*(eye(n)-load_new*load_new')*std_data(i,:)';%Q统计量 end mean_PCA = mean(Q_PCA(1:m)); var_PCA = var(Q_PCA(1:m)); Q_PCA_lim = var_PCA/(2*mean_PCA)*chi2inv(0.99,2*mean_PCA.^2/var_PCA)*ones(length(Q_PCA),1);%求Q统计量的控制限 figure() plot(1:m,Q_PCA,'-r'),title('SPE') xlabel('Sample Number'),ylabel('SPE') hold on plot(Q_PCA_lim(:,1),'b--') hold on plot(1:m,Q_PCA,'og') hold off variance_new = variance(1:r); for i = 1:m T2_PCA(i) = score_new(i,:)*inv(diag(variance_new))*score_new(i,:)'*m;%T2统计量 end T2_PCA_lim = r*(m^2-1)/(m*(m-r))*finv(0.99,r,m-r)*ones(length(T2_PCA),1);%T2统计量的控制限 figure() plot(1:m,T2_PCA,'-r'),title('Hotelling T2') xlabel('Sample Number'),ylabel('T2') hold on plot(T2_PCA_lim(:,1),'b--') hold on plot(1:m,T2_PCA,'og') hold off data_new = score_new*load_new'; figure() plot(std_data(:,1)) hold on plot(data_new(:,1),'r'),title('1st Principal Component') xlabel('Sample Number'),ylabel('Sample Value') legend('original std-data','PCA-SVD data') hold off figure() plot(part2); title('Data for PCA example, Part 2') xlabel('Time') ylabel('Variable Values') [ms,ns] = size(part2); std_data2 = scale(part2,mean_part1,std_part1); score2 = std_data2*load_new; resmat = std_data2'-load_new*score2'; scllim = [1 ms]; res = (sum(resmat.^2))'; q = Q_PCA_lim(1,1); qlim = [q q]; figure() plot(1:ms,res,'-r',scllim,qlim,'--b') hold on plot(1:ms,res,'+g') hold off title('New Sample Residuals with Limits from Old Model') xlabel('Sample Number') ylabel('Residual') tsqvals = sum((score2.^2*inv(diag(ssq(1:r,2))))')'; t = T2_PCA_lim(1,1); tlim = [t t]; figure() plot(1:ms,tsqvals,'-r',scllim,tlim,'--b') hold on plot(1:ms,tsqvals,'+g') hold off title('New Samples Value of T^2 with 99% Limit From Old Model') xlabel('Sample Number') ylabel('Value of T^2')