真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析.docx
《真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析.docx》由会员分享，可在线阅读，更多相关《真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析.docx（18页珍藏版）》请在冰豆网上搜索。
真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析
偏最小二乘回归是一种新型的多元统计数据分析方法，它与1983年由伍德和阿巴诺等人首次提出。
近十年来，它在理论、方法和应用方面都得到了迅速的发展。
密西根大学的弗耐尔教授称偏最小二乘回归为第二代回归分析方法。
偏最小二乘回归方法在统计应用中的重要性主要的有以下几个方面:
（1）偏最小二乘回归是一种多因变量对多自变量的回归建模方法。
（2）偏最小二乘回归可以较好地解决许多以往用普通多元回归无法解决的
问题。
在普通多元线形回归的应用中，我们常受到许多限制。
最典型的问题就是自变量之间的多重相关性。
如果采用普通的最小二乘方法，这种变量多重相关性就会严重危害参数估计，扩大模型误差，并破坏模型的稳定性。
变量多重相关问题十分复杂，长期以来在理论和方法上都未给出满意的答案，这一直困扰着从事
实际系统分析的工作人员。
在偏最小二乘回归中开辟了一种有效的技术途径，它
利用对系统中的数据信息进行分解和筛选的方式，提取对因变量的解释性最强的综合变量，辨识系统中的信息与噪声，从而更好地克服变量多重相关性在系统建模中的不良作用。
（3）偏最小二乘回归之所以被称为第二代回归方法，还由于它可以实现多
种数据分析方法的综合应用0
偏最小二乘回归=多元线性回归分析+典型相关分析+主成分分析
由于偏最小二乘回归在建模的同时实现了数据结构的简化，因此，可以在维平面图上对多维数据的特性进行观察，这使得偏最小二乘回归分析的图形功能十分强大。
在一次偏最小二乘回归分析计算后，不但可以得到多因变量对多自变量的回归模型，而且可以在平面图上直接观察两组变量之间的相关关系，以及观察样本点间的相似性结构。
这种高维数据多个层面的可视见性，可以使数据系统的分析内容更加丰富，同时又可以对所建立的回归模型给予许多更详细深入的实际解释。
偏最小二乘回归的建模策略原理方法
1.1建模原理
设有q个因变量｛yi,...,yq｝和p自变量｛xi,...,xp｝。
为了研究因变量和自变量
的统计关系,我们观测了n个样本点,由此构成了自变量与因变量的数据表
X=｛xi,…,Xp｝和.Y=｛yi,…,yq｝。
偏最小二乘回归分别在X与丫中提取出成分ti和
Ui（也就是说，ti是Xi,...,Xp的线形组合,Ui是yi,…,yq的线形组合）.在提取这
两个成分时,为了回归分析的需要,有下列两个要求:
（1）t1和u1应尽可能大地携带他们各自数据表中的变异信息
（2）t1与u1的相关程度能够达到最大。
这两个要求表明，ti和ui应尽可能好的代表数据表X和丫,同时自变量的成分
ti对因变量的成分ui又有最强的解释能力。
在第一个成分ti和ui被提取后，偏最小二乘回归分别实施X对ti的回归
以及丫对U1的回归。
如果回归方程已经达到满意的精度，贝U算法终止；否则将利用X被ti解释后的残余信息以及丫被ti解释后的残余信息进行第二轮的
成分提取。
如此往复，直到能达到一个较满意的精度为止。
若最终对X共提取
1.2
X经标准化处理后的数
计算方法推导
为了数学推导方便起见,首先将数据做标准化处理。
opnp
据矩阵记为Eo=（Eoi,…，Eop）np，Yj经标准化处理后的数据矩阵记为
F0=（F0i，…，F0q）nP
第一步记ti是Eo的第一个成分，Wi是Eo的第一个轴，它是一个单位向量,
既||w1||=1。
记ui是Fo的第一个成分，u^Foc，。
C1是Fo的第一个轴，并且||ci||=1o
如果要tiUi能分别很好的代表X与丫中的数据变异信息，根据主成分分
析原理，应该有
Var（u1）max
Var（t1）max
另一方面，由于回归建模的需要，又要求ti对ui有很大的解释能力，有典型相关分析的思路，ti与ui的相关度应达到最大值，既
r（t1u1）max
因此，综合起来，在偏最小二乘回归中，我们要求ti与Ui的协方差达到最大，既
Cov（t1
U1）=#Var（tJVar（U1）r（t1，ujmax
正规的数学表述应该是求解下列优化问题，既
max〈E0W1,F0C1
s.tW1W11
C1C11
因此，将在||wi||2=1和||C1||2=1的约束条件下，去求（WE0F0C1）的最大
值。
如果采用拉格朗日算法，记
2（C1C1—1）
s=W1E0FoC1—1（W1W1—1）—
对s分别求关于W1C11和2的偏导并令之为零，有
s=-（c1c1-1）=0
2
由式（1-2）~（1-5），可以推出
2122wlE'oFoqEoW—FoS
记12122w'1E'oFoC1,所以，1正是优化问题的目标函数值.
把式（1-2）和式（1-3）写成
将式（1-7）代入式（1-6）,有
同理，可得
可见,W1是矩阵E'oFoF'oEo的特征向量，对应的特征值为12.1是目标函数值，它要
求取最大值，所以,W1是对应于E'oFoF'oE。
矩阵最大特征值的单位特征向量.而另
II._2、.
一万面，C1是对应于矩阵FoEoEoFo最大特征值1的单位特征向量.
求得轴W1和C1后,即可得到成分
t1E0W1
U1F0C1
然后,分别求Eo和Fo对t1,U1的三个回归方程
Eo
"1
E1
（1-10）
Fo
F1
（1-11）
Fo
t1r'1
F1
（1-12）
式中，回归系数向量是
而E1,F1,F1分别是三个回归方程的残差矩阵.
第二步用残差矩阵E1和F1取代Eo和Fo，然后，求第二个轴W2和C2以及第二个成分t2,U2,有
t2=E1w2
U2=F1C2
2t2,u2w2E1F1C2
C2是对应于矩阵
IIQ
W2是对应于矩阵E1F1F1E1最大特征值2的特征值,
F1t2
||t2||2
因此，有回归方程
由于,t1,,tA均可以表示成Eo1,,Eop的线性组合，因此，式（1-17）还可以还原
FAk是残差距阵Fa的第k列。
1.3交叉有效性
面要讨论的问题是在现有的数据表下,如何确定更好的回归方程。
在许多
情形下,偏最小二乘回归方程并不需要选用全部的成分t1,,tA进行回归建模,而
是可以象在主成分分析一样,采用截尾的方式选择前m个成分
（mA,A秩（X））,仅用这m个后续的成分t1,,tm就可以得到一个预测性较好的模型。
事实上,如果后续的成分已经不能为解释F0提供更有意义的信息时,采用
过多的成分只会破坏对统计趋势的认识,引导错误的预测结论。
在多元回归分析
一章中,我们曾在调整复测定系数的内容中讨论过这一观点。
面的问题是怎样来确定所应提取的成分个数。
在多元回归分析中,曾介绍过用抽样测试法来确定回归模型是否适于预测应用。
我们把手中的数据分成两部分:
第一部分用于建立回归方程,求出回归系数估计量
bB,拟合值?
B以及残差均方和？
B；再用第二部分数据作为实验点，代入刚才所求
得的回归方程，由此求出?
T和?
T。
一般地，若有?
T?
B，则回归方程会有更好的预
测效果。
若?
T?
B,则回归方程不宜用于预测。
在偏最小二乘回归建模中,究竟应该选取多少个成分为宜,这可通过考察增加
一个新的成分后,能否对模型的预测功能有明显的改进来考虑。
采用类似于抽样
测试法的工作方式,把所有n个样本点分成两部分:
第一部分除去某个样本点i的所有样本点集合（共含n-1个样本点）,用这部分样本点并使用h个成分拟合一个回归方程;第二部分是把刚才被排除的样本点i代入前面拟合的回归方程，得到yj在样本点i上的拟合值?
hj（i）。
对于每一个i=1,2,…,n,重复上述测试，则可以定义yj
的预测误差平方和为PRESShj,有
（1-18）
n2
PRESShj（yijy?
hj（i））
i1
定义丫的预测误差平方和为PRESSh，有
（1-19）
P
PRESSPRESSj
j1
显然,如果回归方程的稳健性不好，误差就很大,它对样本点的变动就会十分敏感，这种扰动误差的作用，就会加大PRESSh的值。
另外,再采用所有的样本点，拟合含h个成分的回归方程。
这是，记第i个样本
点的预测值为?
hji，则可以记yj的误差平方和为SShj，有
SShj
（1-20）
定义丫的误差平方和为S&，有
一般说来,总是有PRESSh大于SSh，而SSh则总是小于SSh1。
下面比较SSh1和
PRESSh。
SS.1是用全部样本点拟合的具有h-1个成分的方程的拟合误差
PRESSh增加了一个成分th，但却含有样本点的扰动误差。
如果h个成分的回归方
程的含扰动误差能在一定程度上小于（h-1）个成分回归方程的拟合误差，则认为增加一个成分th,会使预测结果明显提高。
因此我们希望（PRESSh/SSh1）的比值能
越小越好。
在SIMCA-P软件中，指定
（PRESSh/SSh1）0.952
即JPRESS0.95JSS1时,增加成分th就是有益的；或者反过来说,当
JPRESS,0.95JSS7时,就认为增加新的成分th,对减少方程的预测误差无明显
的改善作用.
另有一种等价的定义称为交叉有效性。
对每一个变量yk,定义
对于全部因变量丫,成分th交叉有效性定义为
用交叉有效性测量成分th对预测模型精度的边际贡献有如下两个尺度。
见，q20.0975与（PRESSh/SSh1）0.952是完全等价的决策原则。
可以考虑增加成分th是明显有益的。
明确了偏最小二乘回归方法的基本原理、方法及算法步骤后，我们将做实证分析。
functionw=maxdet（A）%求矩阵的最大特征值
[v,d]=eig（A）;
[n,p]=size（d）;
d1=d*ones（p,1）;
d2=max（d1）;
i=find（d1==d2）;
w=v（:
i）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[c,m,v]=norm1（C）%对数据进行标准化处理
[n,s]=size（C）;
fori=1:
n
forj=1:
s
c（i,j）=（C（i,j）-mean（C（:
j）））/sqrt（cov（C（:
j）））;
endend
m=mean（C）;
forj=1:
s
v（1,j）=sqrt（cov（C（:
j）））;
end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[t,q,w,wh,f0,FF]=fun717（px,py,C）
%px
自变量的输入个数
%py
输入因变量的个数。
%C
输入的自变量和因变量组成的矩阵
%t
提取的主成分
%q
为回归系数。
%w
最大特征值所对应的特征向量。
E0=c（:
1:
px）;
F0=c（:
px+1:
px+py）;
A=E0'*F0*F0'*E0;
提取主成分
t（:
1）=E0*w（:
1）;
E（:
1:
px）=E0-t（:
1）*（E0'*t（:
1）/（t（:
1）'*t（:
1）））';
获得回归系数
p（:
1:
px）=（E0'*t（:
1）/（t（:
1）'*t（:
1）））';
fori=0:
px-2
B（:
px*i+1:
px*i+px）=E（:
px*i+1:
px*i+px）'*F0*F0'*E（:
px*i+1:
px*i+px）
w（:
i+2）=maxdet（B（:
px*i+1:
px*i+px））;
%maxdet为求最大特征值的函数
t（:
i+2）=E（:
px*i+1:
px*i+px）*w（:
i+2）;
p（:
px*i+px+1:
px*i+2*px）=（E（:
px*i+1:
px*i+px）'*t（:
i+2）/（t（:
i+2）'*t（:
i+2）））';
E（:
px*i+px+1:
px*i+2*px）=E（:
px*i+1:
px*i+px）-t（:
i+2）*（E（:
px*i+1:
px*i+px）'*t（:
i+2）/（t（:
i+2）'*t（:
i+2）））';
end
fors=1:
px
求回归系数
%no
q（:
s）=p（1,px*（s-1）+1:
px*s）';
end
[n,d]=size（q）;
forh=1:
px
iw=eye（d）;
forj=1:
h-1
iw=iw*（eye（d）-w（:
j）*q（:
j）'）;
endwh（:
h）=iw*w（:
h）;
end
forj=1:
py
zr（j,:
）=（regress1（y（:
j）,t））';
end
forj=1:
px
fori=1:
py
生成标准化变量的方程的系数矩阵
w1=wh（:
1:
j）;
zr1=（zr（i,1:
j））';
f0（i,:
j）=（w1*zr1）';
end
[normxy,meanxy,covxy]=norm1（C）;
rmxy标准化后的数据矩阵
%meanx每一列的均值
%covxy每一列的方差
ccxx=ones（py,1）*meanxy（1,1:
px）;
ccy=（covxy（1,px+1:
px+py））'*ones（1,px）;
ccx=ones（py,1）*（covxy（1,1:
px））;
ff=ccy.*f0（:
:
j）./ccx;
生成
fff=-（sum（（ccy.*ccxx.*f0（:
:
j）./ccx）'）-meanxy（1,px+1:
px+py））';
FF（:
:
j）=[fff,ff];
原始变量方程的常数项和系数矩阵
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[r,Rdyt,RdYt,RdYtt,Rdytt,VIP]=fun8y（px,py,c）
X=c（:
1:
px）;
Y=c（:
px+1:
px+py）;
x=norm1（X）;
y=norm1（Y）;
[t,q,w]=fun717（px,py,[X,Y]）;
r1=corrcoef（[y,t]）;
r=r1（py+1:
px+py,1:
py）';
Rdyt=r.A2;
RdYt=mean（Rdyt）form=1:
px
RdYtt（1,m）=sum（RdYt（1,1:
m）'）;
end
form=1:
py
Rdytt（j,m）=sum（Rdyt（j,1:
m）'）;
end
end
forj=1:
px
form=1:
px
Rd（j,m）=RdYt（1,1:
m）*（（w（j,1:
m）A2）'）;
end
endforj=1:
px
VIP（j,:
）=sqrt（（px*ones（1,px）./RdYtt）.*Rd（j,:
））;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[r,Rdxt,RdXt,RdXtt,Rdxtt]=fun8x（px,py,c）
X=c（:
1:
px）;
Y=c（:
px+1:
px+py）;
x=norm1（X）;
y=norm1（Y）;
[t,q,w]=fun717（px,py,[X,Y]）;
r1=corrcoef（[x,t]）;
r=r1（px+1:
px+px,1:
px）';
Rdxt=r.A2;
RdXt=mean（Rdxt）;
form=1:
px
RdXtt（1,m）=sum（RdXt（1,1:
m）'）;
endforj=1:
px
Rdxtt（j,m）=sum（Rdxt（j,1:
m）'）;
end
end
%forj=1:
px
form=1:
px
Rd（j,m）=RdXt（1,1:
m）*（（w（j,1:
m）A2）'）;
end
endforj=1:
px
VIP（j,:
）=sqrt（（px*ones（1,px）./RdYtt）.*Rd（j,:
））;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[t,u]=TU（px,py,C）
%t提取的自变量的主成分
%u提取的因变量的主成分
c=norm1（C）;
y=c（:
px+1:
px+py）;
E0=c（:
1:
px）;
F0=c（:
px+1:
px+py）;
A=E0'*F0*F0'*E0;
w（:
1）=maxdet（A）;
t（:
1）=E0*w（:
1）;
B=F0'*E0*E0'*F0;
cc（:
1）=maxdet（B）;
u（:
1）=F0*cc（:
1）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functiondrew（px,py,c）
X=c（:
1:
px）;
Y=c（:
px+1:
px+py）;
[line,l]=size（Y）;
[t,q,w,wh,f0,FF]=fun717（px,py,c）;
YY=X*FF（:
2:
px+1,3）'+ones（line,1）*FF（:
1,3）';
subplot（1,1,1,1）
bar（f0（:
:
3））
legend（'SG','TZBFB','FHL','JK','HPZD','JPZD','TZ','ZG','GPK'）
gridon
plot（YY（:
4）,Y（:
4）,'+'）;
lsline
fori=1:
py
v=mod（i,4）;
d=（i-v）/4;
subplot（2,2,v,d+1）plot（YY（:
i）,Y（:
i）,'*'）;
lsline
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[Qhj,Qh,prey]=crossval7（px,py,c）
%px是自变量的个数；
%py是因量
PRESShj=zeros（px,py）;
X=c（:
1:
px）;
Y=c（:
px+1:
px+py）;
x=norm1（X）;
y=norm1（Y）;
[line,row]=size（x）;
forj=1:
line
newx=X;
newy=Y;
newx（j,:
）=[];
newy（j,:
）=[];
[t,p0,w,wh,f0,FF]=fun717（px,py,[newx,newy]）;
prey（j,:
h）=X（j,:
）*FF（:
2:
px+1,h）'+FF（:
1,h）';
end
PRESShj（h,:
）=sum（（Y-prey（:
，:
6）八2）;
end
PRESSh=（sum（PRESShj'））';
forh=1:
px
[t1,p0,w,wh,f0,FF]=fun717（px,py,c）;
prey2（:
:
h）=X（:
:
）*FF（:
2:
px+1,h）'+ones（line,1）*F
F（:
1,h）';
SShj（h,:
）=sum（（Y-prey2（:
:
h））A2）;
end
SSh=（sum（SShj'））';
Qhj=ones（px-1,py）-PRESShj（2:
px,:
）./SShj（1:
px-1,:
）;%错位
Qh=ones（px-1,1）-PRESSh（2:
px,1）./SSh（1:
px-1,1）;
升级会员 