软测量作业matlab_用偏最小二乘法PLS求解回归方程

题目：

有7个自变量 ，因变量记为 ，即：

——直接蒸馏成分； ——聚合物；
——重整汽油； ——烷基化物；
——原油热裂化油； ——天然香精；
——原油催化裂化油； ——原辛烷值；
表1给出了12种混合物关于这8个变量的观测数据，请利用PLS建立 对 的回归方程，以确定7种构成元素 对 的影响。

表1 观测数据表
NO x_1 x_2 x_3 x_4 x_5 x_6 x_7 y
1 0.00 0.23 0.00 0.00 0.00 0.74 0.03 98.7
2 0.00 0.10 0.00 0.00 0.12 0.74 0.04 97.8
3 0.00 0.00 0.00 0.10 0.12 0.74 0.04 96.6
4 0.00 0.49 0.00 0.00 0.12 0.37 0.02 92.0
5 0.00 0.00 0.00 0.62 0.12 0.18 0.08 86.6
6 0.00 0.62 0.00 0.00 0.00 0.37 0.01 91.2
7 0.17 0.27 0.10 0.38 0.00 0.00 0.08 81.9
8 0.17 0.19 0.10 0.38 0.02 0.06 0.08 83.1
9 0.17 0.21 0.10 0.38 0.00 0.06 0.08 82.4
10 0.17 0.15 0.10 0.38 0.02 0.10 0.08 83.2
11 0.21 0.36 0.12 0.25 0.00 0.00 0.06 81.4
12 0.00 0.00 0.00 0.55 0.00 0.37 0.08 88.1

代码：

xy=[0.00	0.23	0.00	0.00	0.00	0.74	0.03	98.7
0.00	0.10	0.00	0.00	0.12	0.74	0.04	97.8
0.00	0.00	0.00	0.10	0.12	0.74	0.04	96.6
0.00	0.49	0.00	0.00	0.12	0.37	0.02	92.0
0.00	0.00	0.00	0.62	0.12	0.18	0.08	86.6
0.00	0.62	0.00	0.00	0.00	0.37	0.01	91.2
0.17	0.27	0.10	0.38	0.00	0.00	0.08	81.9
0.17	0.19	0.10	0.38	0.02	0.06	0.08	83.1
0.17	0.21	0.10	0.38	0.00	0.06	0.08	82.4
0.17	0.15	0.10	0.38	0.02	0.10	0.08	83.2
0.21	0.36	0.12	0.25	0.00	0.00	0.06	81.4
0.00	0.00	0.00	0.55	0.00	0.37	0.08	88.1
];
%数据标准化
x0=xy(:,1:end-1);
y0=xy(:,end);
std_xy=std(xy);
mean_xy=mean(xy);
xy_bz=zscore(xy);
e0=xy_bz(:,1:end-1);
f0=xy_bz(:,end); 

[num_row,num_col]=size(xy);
n=num_col-1;
I=eye(n);   

%设置循环，直到满足交叉有效性的条件。
for i=1:n
    X=e0'*f0*f0'*e0;
    [V,D]=eig(X);   %特征向量V，特征值D
    D=diag(D);    %提出特征值
    [max_D,index] = max(D); %提取出最大特征值和其索引
    w(:,i)=V(:,index);    %提出最大特征值对应的特征向量
    w_star(:,i)=I*w(:,i); %创建w*
    
    t(:,i)=e0*w(:,i); %ti的表达式
    alpha=e0'*t(:,i)/(t(:,i)'*t(:,i)); %P1，用来求W*
    I=I*(eye(n)-w(:,i)*alpha'); %计算w 到w*的变换矩阵
    e=e0-t(:,i)*alpha'; %计算残差矩阵
    e0=e;%把新的E1赋值过去迭代计算
    
    %预测误差平方和PRESS
    for j=1:num_row
        t1=t(:,1:i);f1=f0;
        she_t=t1(j,:);she_f=f1(j,:); %把舍去的第j个样本点保存起来
        t1(j,:)=[];f1(j,:)=[]; %删除第j个观测值
        beta1=[t1,ones(num_row-1,1)]\f1; %求回归分析的系数
        beta1(end,:)=[]; %删除回归分析的常数项
        cancha=she_f-she_t*beta1; %求残差向量
        p_i(j)=sum(cancha.^2);   %误差平方和
    end
    p(i)=sum(p_i);
    
    %误差平方和 SS(h-1)
    beta=[t(:,1:i),ones(num_row,1)]\f0; %求回归方程的系数
    beta(end,:)=[]; %删除回归分析的常数项
    cancha=f0-t(:,1:i)*beta; %预测误差
    ss(i)=sum(sum(cancha.^2)); %预测误差平方和 
    
    % 交叉有效性检验,EFF
    if i>1
        EFF(i)=1-p(i)/ss(i-1);  % Q_h2为1-PRESS(h)/SS(h-1)
    else
        EFF(1)=1;
    end
    if EFF(i)<0.0975
        fprintf('提出了%d个成分',i);
        EFF(i)
        r=i;%用r记录i
        break
    end
end

% 求原始回归方程的系数
beta_bz=[t(:,1:r),ones(num_row,1)]\f0;%根据提出的成分建立标准后的回归方程
beta_bz(end,:)=[];%去掉多余的
xishu=w_star(:,1:r)*beta_bz; %F=AB+FA

C=mean_xy(end)-mean_xy(1:n)./std_xy(1:n)*std_xy(end)*xishu(:,1)    %原始数据的常数项
beta_x(:,1)=xishu(:,1)./std_xy(1:n)'*std_xy(end) %原始数据的系数

%对比真实输出与模型估计输出,计算两组数据的残差平方和RSS和回归平方和ESS。
Y=xy(:,end);
y_mean=mean(Y);
y_hat=xy(:,1:end-1)*beta_x+C
RSS=0;
ESS=0;
TSS=0;
for i=1:12
   RSS=RSS+(Y(i)-y_hat(i))^2 ;
   ESS=ESS+(y_hat(i)-y_mean)^2;
end
TSS=ESS+RSS;
r2=ESS/TSS

%绘制图像
plot(80:100,80:100,y_hat,Y,'*')