真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析Word格式.docx
《真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析Word格式.docx》由会员分享,可在线阅读,更多相关《真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析Word格式.docx(18页珍藏版)》请在冰点文库上搜索。
为了数学推导方便起见,首先将数据做标准化处理。
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'
oE。
矩阵最大特征值的单位特征向量.而另
II._2、.
一万面,C1是对应于矩阵FoEoEoFo最大特征值1的单位特征向量.
求得轴W1和C1后,即可得到成分
t1E0W1
U1F0C1
然后,分别求Eo和Fo对t1,U1的三个回归方程
Eo
"
1
E1
(1-10)
Fo
F1
(1-11)
t1r'
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,则回归方程会有更好的预
测效果。
若?
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);
v(1,j)=sqrt(cov(C(:
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(:
px)=E0-t(:
1)*(E0'
*t(:
1)/(t(:
1)'
1)))'
;
获得回归系数
p(:
px)=(E0'
fori=0:
px-2
B(:
px*i+1:
px*i+px)=E(:
px*i+px)'
*E(:
px*i+px)
w(:
i+2)=maxdet(B(:
px*i+px));
%maxdet为求最大特征值的函数
i+2)=E(:
px*i+px)*w(:
i+2);
px*i+px+1:
px*i+2*px)=(E(:
i+2)/(t(:
i+2)'
*t(:
i+2)))'
px*i+2*px)=E(:
px*i+px)-t(:
i+2)*(E(:
px*i+px)'
end
fors=1:
px
求回归系数
%no
q(:
s)=p(1,px*(s-1)+1:
px*s)'
[n,d]=size(q);
forh=1:
iw=eye(d);
h-1
iw=iw*(eye(d)-w(:
j)*q(:
j)'
);
endwh(:
h)=iw*w(:
h);
py
zr(j,:
)=(regress1(y(:
j),t))'
fori=1:
生成标准化变量的方程的系数矩阵
w1=wh(:
j);
zr1=(zr(i,1:
j))'
f0(i,:
j)=(w1*zr1)'
[normxy,meanxy,covxy]=norm1(C);
rmxy标准化后的数据矩阵
%meanx每一列的均值
%covxy每一列的方差
ccxx=ones(py,1)*meanxy(1,1:
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:
FF(:
j)=[fff,ff];
原始变量方程的常数项和系数矩阵
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[r,Rdyt,RdYt,RdYtt,Rdytt,VIP]=fun8y(px,py,c)
X=c(:
Y=c(:
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:
RdYtt(1,m)=sum(RdYt(1,1:
m)'
form=1:
Rdytt(j,m)=sum(Rdyt(j,1:
Rd(j,m)=RdYt(1,1:
m)*((w(j,1:
m)A2)'
endforj=1:
VIP(j,:
)=sqrt((px*ones(1,px)./RdYtt).*Rd(j,:
));
function[r,Rdxt,RdXt,RdXtt,Rdxtt]=fun8x(px,py,c)
r1=corrcoef([x,t]);
r=r1(px+1:
px+px,1:
px)'
Rdxt=r.A2;
RdXt=mean(Rdxt);
RdXtt(1,m)=sum(RdXt(1,1:
Rdxtt(j,m)=sum(Rdxt(j,1:
%forj=1:
Rd(j,m)=RdXt(1,1:
function[t,u]=TU(px,py,C)
%t提取的自变量的主成分
%u提取的因变量的主成分
c=norm1(C);
y=c(:
1)=maxdet(A);
B=F0'
*E0*E0'
*F0;
cc(:
1)=maxdet(B);
u(:
1)=F0*cc(:
functiondrew(px,py,c)
[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
v=mod(i,4);
d=(i-v)/4;
subplot(2,2,v,d+1)plot(YY(:
i),Y(:
i),'
*'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[Qhj,Qh,prey]=crossval7(px,py,c)
%px是自变量的个数;
%py是因量
PRESShj=zeros(px,py);
[line,row]=size(x);
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(:
px+1,h)'
+FF(:
1,h)'
PRESShj(h,:
)=sum((Y-prey(:
,:
6)八2);
PRESSh=(sum(PRESShj'
))'
[t1,p0,w,wh,f0,FF]=fun717(px,py,c);
prey2(:
h)=X(:
+ones(line,1)*F
F(:
SShj(h,:
)=sum((Y-prey2(:
h))A2);
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);