matlab进行pca分析,matlab 进行PCA分析

% This program is used to explain how PCA works and the idea behind the

% method. This example is given in the ISA book (Chapter 8, simple example)

% The process is composed of a valve a thermometer and a fermenter. The

% valve is used to control the flow rate of cooling water. The thermometer

% is used to measure the temperature of the fermenter.

% The steady state of the system is: valve=10% max. Temperature T=25 D.

% YANG ZHANG Email: yz@che.utexas.edu

% 2006-April, UTexas-Austin

clear all

% Give valve set points and fermenter temperature.

% Valve=10+rand(1,10)-0.5; % ten value of Valve around mean.

% Temp=25-(Valve-10)+(rand(1,10)-0.5)*0.5; % 10 value of Temp: T=2.5*valve + noise

Temp = [2.5000e+001 2.5000e+001 2.5000e+001 2.5300e+001 2.5200e+001 2.5200e+001 2.4800e+001 2.4700e+001 2.5100e+001 2.4700e+001];

Valve = [1.0100e+001 9.9000e+000 1.0000e+001 9.8000e+000 9.9000e+000 9.8000e+000 1.0100e+001 1.0300e+001 1.0000e+001 1.0100e+001];

% Define database

x=[Valve;Temp]';

figure(1)

plot(x(:,1),x(:,2),'ko')

xlabel('Valve Position (%)')

ylabel('Fermenter Temperature (C)')

% Calculate the covariance of dataset x (relationship between Valve and

% Temperature)

covariance=cov(x);

% Singular Value Decomposition x=T*E*P

% E is the eignvalue of x; U is orthnormal; V=U'.

%[U E V]=svd(covariance);

%U=U'; V=V';

[U E] = eig(covariance);

E = diag([E(2,2),E(1,1)]); % arrange the Eignvalue in a decending order

U = [U(:,2),U(:,1)]; % arrange the Eignvector to make it compatable with Eignvalue.

% Mean center the raw data

Valvebar=x(:,1)-mean(x(:,1));

Tempbar=x(:,2)-mean(x(:,2));

xbar=[Valvebar,Tempbar];

% Principal axis rotation of the covariance matrix.

z=U'*xbar';

z=z';

%covz is the covariance matrix of principal components which is equal to E

%matrix. The variance of transformed variable (z1 & z2) will have variance

%0.2283 and 0.0123 respectively.

covz=U'*covariance*U;

% The first column of U factor is - and + with nearly equal value which

% means the first principal component is related to variability which both

% measurements have difference. The second column is both positive which

% means the 2nd PC is concentrate on the varaiability of the common between

% the two.

% scaling of PCS.

for i=1:2

Vs(:,i)=sqrt(E(i,i))*U(:,i);

Ws(:,i)=U(:,i)/sqrt(E(i,i));

end

% after rescale: Vs'*Vs=E; V'*COV*Vs=E^2; Ws'*Ws=inv(E); Ws'*COV*Ws=I

y=Ws'*xbar';

%define T2

T2=diag(y'*y);

% find the control limit of T2

p=2; %PC number

n=10; % sample number

T2Upper=p*(n-1)/(n-p)*finv(0.95,p,n-p);

% draw the T2 control charts

figure(2)

plot(T2,'ko')

hold on

plot(T2Upper*ones(1,13),'k')

xlabel('Sample Number')

ylabel('T2 Score')

hold off

% Draw the control ellipse

% s1^2*s2^2/(s1^2*s2^2-s12^2)*[(x1-x1mean)^2/s1^2+(x2-x2mean)^2/s2^2-2*s12(

% x1-x1mean)(x2-x2mean)/s1^2/s2^2]=T2^2(upper control limit)

% S=covariance;

% a=S(1)*S(4)/(S(1)*S(4)-S(2)^2);

% b=xbar(:,1).^2/S(1)+xbar(:,2).^2/S(4)-2*S(2)*xbar(:,1).*xbar(:,2)./S(1)/S(4);

% error detection

% two new points is added to be detected which are [10.1,25.2];[9.9;24];

% plot two values with the original 10.

Valve2=[10.2,9.9];Temp2=[25.1,25.5];

figure(1)

hold on

plot(Valve2,Temp2,'k*')

hold off

% plot the control charts for valve and temperature separately.

figure(3)

hold on

plot(1:10,Valve,'ko')

plot(11:12,Valve2,'k*')

plot(0:13, mean(Valve)*ones(1,14),'k')

plot(0:13, (mean(Valve)-std(Valve)*2)*ones(1,14),'k')

plot(0:13, (mean(Valve)+std(Valve)*2)*ones(1,14),'k')

xlabel('Sample Number')

ylabel('Valve Position (%)')

hold off

figure(4)

hold on

plot(1:10,Temp,'ko')

plot(11:12,Temp2,'k*')

plot(0:13, mean(Temp)*ones(1,14),'k')

plot(0:13, (mean(Temp)+std(Temp)*2)*ones(1,14),'k')

plot(0:13, (mean(Temp)-std(Temp)*2)*ones(1,14),'k')

xlabel('Sample Number')

ylabel('Temperature (C)')

hold off

% transform new observations to PC axis

y2=Ws'*([Valve2-mean(Valve);Temp2-mean(Temp)]);

% plot the control charts for principal components

figure(5)

hold on

plot(1:10,y(1,:),'ko')

plot(11:12,y2(1,:),'k*')

plot(0:13, mean(y(1,:))*ones(1,14),'k')

plot(0:13, -2*ones(1,14),'k')

plot(0:13, 2*ones(1,14),'k')

xlabel ('Sample Number')

ylabel ('Principal Component 1')

hold off

figure(6)

hold on

plot(1:10,y(2,:),'ko')

plot(11:12,y2(2,:),'k*')

plot(0:13,mean(y(2,:))*ones(1,14),'k')

plot(0:13, -2*ones(1,14),'k')

plot(0:13, 2*ones(1,14),'k')

xlabel ('Sample Number')

ylabel ('Principal Component 2')

hold off

% plot the T2 chart for the two new points

figure(2)

T22=diag(y2'*y2);

hold on

plot([11,12],T22,'k*')

hold off

% draw the control ellipse

ra = sqrt(T2Upper*E(1,1));

rb = sqrt(T2Upper*E(2,2));

V=[sqrt(E(1,1))*U(:,1),sqrt(E(2,2))*-U(:,2)];

k1 = V(2,1)/V(1,1);

k2 = V(2,2)/V(1,2);

ang = atan(k1); %atan(k1);

x0 = mean(x(:,1));

y0 = mean(x(:,2));

% i = 9.5:0.01:10.5;

% s = 9.72:0.01:10.3;

% j = k1*i + 42.6303;

% l = k2*s + 19.3279;

figure (7)

h=ellipse(ra,rb,ang,x0,y0,'k',1000);

hold on

plot(x(:,1),x(:,2),'ko')

plot(Valve2, Temp2,'k*')

% plot(i,j)

% plot(s,l)

xlabel ('Valve Position')

ylabel ('Temperature')

hold off