第三十章偏最小二乘回归 在实际问题中,经常遇到需要研究两组多重相关变量间的相互依赖关系,并研究用 组变量(常称为自变量或预测变量)去预测另一组变量(常称为因变量或响应变量), 除了最小二乘准则下的经典多元线性回归分析(MLR),提取自变量组主成分的主成 分回归分析(PCR)等方法外,还有近年发展起来的偏最小二乘(PLS)回归方法 偏最小二乘回归提供一种多对多线性回归建模的方法,特别当两组变量的个数很 多,且都存在多重相关性,而观测数据的数量(样本量)又较少时,用偏最小二乘回归 建立的模型具有传统的经典回归分析等方法所没有的优点 偏最小二乘回归分析在建模过程中集中了主成分分析,典型相关分析和线性回归分 析方法的特点,因此在分析结果中,除了可以提供一个更为合理的回归模型外,还可以 同时完成一些类似于主成分分析和典型相关分析的研究内容,提供更丰富、深入的一些 信息。 本章介绍偏最小二乘回归分析的建模方法;通过例子从预测角度对所建立的回归模 型进行比较。 §1偏最小二乘回归 考虑P个变量y,y2…,yn与m个自变量x,x2…,xm的建模问题。偏最小二乘 回归的基本作法是首先在自变量集中提出第一成分1(1是x1,…,xm的线性组合,且 尽可能多地提取原自变量集中的变异信息);同时在因变量集中也提取第一成分l1 并要求h与相关程度达到最大。然后建立因变量y,…,y与4的回归,如果回归方 程已达到满意的精度,则算法中止。否则继续第二对成分的提取,直到能达到满意的精 度为止。若最终对自变量集提取厂个成分1,12…1,偏最小二乘回归将通过建立 n,…y与4,42…的回归式,然后再表示为y…,y与原自变量的回归方程式 即偏最小二乘回归方程式 为了方便起见,不妨假定P个因变量y,…,yn与m个自变量x,…,xm均为标准 化变量。因变量组和自变量组的n次标准化观测数据阵分别记为 y F Eo 偏最小二乘回归分析建模的具体步骤如下:
-531- 第三十章 偏最小二乘回归 在实际问题中，经常遇到需要研究两组多重相关变量间的相互依赖关系，并研究用 一组变量（常称为自变量或预测变量）去预测另一组变量（常称为因变量或响应变量）， 除了最小二乘准则下的经典多元线性回归分析（MLR），提取自变量组主成分的主成 分回归分析（PCR）等方法外，还有近年发展起来的偏最小二乘（PLS）回归方法。 偏最小二乘回归提供一种多对多线性回归建模的方法，特别当两组变量的个数很 多，且都存在多重相关性，而观测数据的数量（样本量）又较少时，用偏最小二乘回归 建立的模型具有传统的经典回归分析等方法所没有的优点。 偏最小二乘回归分析在建模过程中集中了主成分分析，典型相关分析和线性回归分 析方法的特点，因此在分析结果中，除了可以提供一个更为合理的回归模型外，还可以 同时完成一些类似于主成分分析和典型相关分析的研究内容，提供更丰富、深入的一些 信息。 本章介绍偏最小二乘回归分析的建模方法；通过例子从预测角度对所建立的回归模 型进行比较。 §1 偏最小二乘回归 考虑 p 个变量 p y , y , , y 1 2 " 与 m 个自变量 m x , x , , x 1 2 " 的建模问题。偏最小二乘 回归的基本作法是首先在自变量集中提出第一成分 1t （ 1t 是 m x , , x 1 " 的线性组合，且 尽可能多地提取原自变量集中的变异信息）；同时在因变量集中也提取第一成分 1 u ， 并要求 1t 与 1 u 相关程度达到最大。然后建立因变量 p y , , y 1 " 与 1t 的回归，如果回归方 程已达到满意的精度，则算法中止。否则继续第二对成分的提取，直到能达到满意的精 度为止。若最终对自变量集提取 r 个成分 r t ,t , ,t 1 2 " ，偏最小二乘回归将通过建立 p y , , y 1 " 与 r t ,t , ,t 1 2 " 的回归式，然后再表示为 p y , , y 1 " 与原自变量的回归方程式， 即偏最小二乘回归方程式。 为了方便起见，不妨假定 p 个因变量 p y , , y 1 " 与 m 个自变量 m x , , x 1 " 均为标准 化变量。因变量组和自变量组的n 次标准化观测数据阵分别记为 ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = n np p y y y y F " # # " 1 11 1 0 ， ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = n nm m x x x x E " # # " 1 11 1 0 偏最小二乘回归分析建模的具体步骤如下：
(1)分别提取两变量组的第一对成分,并使之相关性达最大 假设从两组变量分别提出第一对成分为1和l1,1是自变量集X=(x1,…,xn)的 线性组合:1=W1x+…+1nxn=w1X,1是因变量集Y=(y,…,yn)的线性组 合:4=1y+…+Vpy=v1Y。为了回归分析的需要,要求 ①t1和1各自尽可能多地提取所在变量组的变异信息 ②1和1的相关程度达到最大。 由两组变量集的标准化观测数据阵E0和F0,可以计算第一对成分的得分向量,记 为1和 xInWu E Fov 第一对成分1和1的协方差Cov(1,1)可用第一对成分的得分向量1和的内积 来计算。故而以上两个要求可化为数学上的条件极值问题: ==WEFx→max =1,rn=|n 利用 Lagrange乘数法,问题化为求单位向量w1和v1,使B1=wEFV1→最大。问 题的求解只须通过计算m×m矩阵M=E6F0FE的特征值和特征向量,且M的最大特 征值为,相应的单位特征向量就是所求的解,而v可由w计算得到n=FE (2)建立y,…,yn对1的回归及x1…,xm对t1的回归。 假定回归模型为
-532- （1）分别提取两变量组的第一对成分，并使之相关性达最大。 假设从两组变量分别提出第一对成分为 1t 和 1 u ，1t 是自变量集 T m X (x , , x ) = 1 " 的 线性组合：t w x w x w XT 1 = 11 1 +"+ 1m m = 1 ， 1 u 是因变量集 T p Y ( y , , y ) = 1 " 的线性组 合：u v y v y v Y T 1 = 11 1 +"+ 1p p = 1 。为了回归分析的需要，要求： ① 1t 和 1 u 各自尽可能多地提取所在变量组的变异信息； ② 1t 和 1 u 的相关程度达到最大。 由两组变量集的标准化观测数据阵 E0 和 F0 ，可以计算第一对成分的得分向量，记 为 1t ˆ 和 1 uˆ ： ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = = 1 11 1 11 1 11 1 1 0 1 ˆ n nm m n m t t w w x x x x t E w # # " # # " ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = = 1 11 1 11 1 11 1 1 0 1 ˆ n np p n p u u v v y y y y u F v # # " # # " 第一对成分 1t 和 1 u 的协方差Cov( , ) 1 1 t u 可用第一对成分的得分向量 1t ˆ 和 1 uˆ 的内积 来计算。故而以上两个要求可化为数学上的条件极值问题： ⎪⎩ ⎪ ⎨ ⎧ = = = = == ⇒ 1, 1 ˆ , ˆ , max 2 1 1 1 2 1 1 1 1 0 1 0 1 1 0 0 1 w w w v v v t u E w Y v w E F x T T T T 利用Lagrange乘数法，问题化为求单位向量 w1和 1 v ，使θ1 = w1 T E0 T F0v1 ⇒ 最大。问 题的求解只须通过计算 m× m 矩阵 M E0 F0F0 E0 T T = 的特征值和特征向量，且 M 的最大特 征值为 2 θ1 ，相应的单位特征向量就是所求的解 w1，而 1 v 可由 w1计算得到 0 0 1 1 1 1 v F E w T θ = 。 （2）建立 p y , , y 1 " 对 1t 的回归及 m x , , x 1 " 对 1t 的回归。 假定回归模型为
alte Fo=uB:+F 其中a1=(a1…a1n),B1=(1…,月n)分别是多对一的回归模型中的参数向量, E1和F1是残差阵。回归系数向量a1,B1的最小二乘估计为 a=E2/ [B=Foi,/l 称a1,B1为模型效应负荷量 (3)用残差阵E和F代替E和F重复以上步骤 记E0=1a1,F0=iB,则残差阵E1=E6-E0,F=F-F0。如果残差阵F 中元素的绝对值近似为0,则认为用第一个成分建立的回归式精度已满足需要了,可以 停止抽取成分。否则用残差阵E1和F代替E和F重复以上步骤即得 W2=(21…,w2m);v2=(n2,…,2n)分别为第二对成分的权数。而 2=E12,i2=Fv2为第二对成分的得分向量。 分别为X,Y的第二对成分的负荷量。这时有 E0=1a1 E F=iB+2月2+F2 (4)设n×m数据阵E0的秩为r≤min(n-1,m),则存在r个成分1,12… 使得 t a+e F=i1B1+…+iB+F 把l=Wk1x1+…+mxn(k=1,2…r),代入Y=1月1+…+1,B,即得P个 因变量的偏最小二乘回归方程式
-533- ⎪⎩ ⎪ ⎨ ⎧ = + = + 0 1 1 1 0 1 1 1 ˆ ˆ F u F E t E T T β α 其中 T m ( , , ) α1 = α11 " α1 ， T p ( , , ) β1 = β11 " β1 分别是多对一的回归模型中的参数向量， E1和 F1是残差阵。回归系数向量 1 1 α , β 的最小二乘估计为 ⎪ ⎩ ⎪ ⎨ ⎧ = = 2 1 0 1 1 2 1 0 1 1 ˆ ˆ ˆ ˆ F t t E t t T T β α ， 称 1 1 α , β 为模型效应负荷量。 （3）用残差阵 E1和 F1代替 E0 和 F0 重复以上步骤。 记 T E t 0 1 1 ˆ = ˆα ， T F t 0 1 1 ˆ = ˆ β ，则残差阵 1 0 0 E = E − Eˆ ， 1 0 0 F = F − Fˆ 。如果残差阵 F1 中元素的绝对值近似为 0，则认为用第一个成分建立的回归式精度已满足需要了，可以 停止抽取成分。否则用残差阵 E1和 F1代替 E0 和 F0 重复以上步骤即得： T w w w m ( , , ) 2 = 21 " 2 ； T p v (v , , v ) 2 = 21 " 2 分别为第二对成分的权数。而 2 1 2 t ˆ = E w ， 2 1 2 uˆ = F v 为第二对成分的得分向量。 2 2 1 2 2 ˆ ˆ t E t T α = ， 2 2 1 2 2 ˆ ˆ t F t T β = 分别为 X ,Y 的第二对成分的负荷量。这时有 ⎪⎩ ⎪ ⎨ ⎧ = + + = + + 0 1 1 2 2 2 0 1 1 2 2 2 ˆ ˆ ˆ ˆ F t t F E t t E T T T T β β α α （4）设 n × m 数据阵 E0 的秩为 r ≤ min(n −1,m) ，则存在 r 个成分 r t ,t , ,t 1 2 " ， 使得 ⎪⎩ ⎪ ⎨ ⎧ = + + + = + + + r T r r T r T r r T F t t F E t t E β β α α ˆ ˆ ˆ ˆ 0 1 1 0 1 1 " " 把 k k km m t = w x +"+ w x 1 1 （ k = 1,2,",r ），代入 r r Y = t1β1 +"+ t β ，即得 p 个 因变量的偏最小二乘回归方程式
y (5)交叉有效性检验。 一般情况下,偏最小二乘法并不需要选用存在的r个成分12…来建立回归 式,而像主成分分析一样,只选用前l个成分(l≤r),即可得到预测能力较好的回归 模型。对于建模所需提取的主成分个数l,可以通过交叉有效性检验来确定 每次舍去第i个观测(i=1,2,…,n),用余下的n-1个观测值按偏最小二乘回归 方法建模,并考虑抽取h个成分后拟合的回归式,然后把舍去的第i个观测点代入所拟 合的回归方程式,得到y(j=1,2…,P)在第i个观测点上的预测值yy(h)。对 i=1,2,…,n重复以上的验证,即得抽取h个成分时第j个因变量y,(=1,2,…,p)的 预测误差平方和为 PRESS/(h)=∑(y-jn/(h)(j=12…P) =(y1…,yn)的预测误差平方和为 PRESSO(h)=∑ PRESS/(h) 另外,再采用所有的样本点,拟合含h个成分的回归方程。这时,记第i个样本点 的预测值为y(h),则可以定义y的误差平方和为 S8(b)=∑(yn-j(h)2 定义Y的误差平方和为 SSh)=∑SS(h) 当 PRESS(h)达到最小值时,对应的h即为所求的成分个数。通常,总有 PRESS(h)大于SS(h),而SS(h)则小于SS(h-1)。因此,在提取成分时,总希望比 值 PRESS(h)/SS(h-1)越小越好;一般可设定限制值为05,即当 534
-534- j j jm m y = a x +"+ a x 1 1 ，（ j = 1,2,",m） （5）交叉有效性检验。 一般情况下，偏最小二乘法并不需要选用存在的 r 个成分 r t ,t , ,t 1 2 " 来建立回归 式，而像主成分分析一样，只选用前l 个成分（l ≤ r ），即可得到预测能力较好的回归 模型。对于建模所需提取的主成分个数l ，可以通过交叉有效性检验来确定。 每次舍去第i 个观测（i = 1,2,", n ），用余下的 n −1个观测值按偏最小二乘回归 方法建模，并考虑抽取h 个成分后拟合的回归式，然后把舍去的第i 个观测点代入所拟 合的回归方程式，得到 y ( j 1,2, , p) j = " 在第i 个观测点上的预测值 ˆ ( ) y(i) j h 。对 i = 1,2,", n 重复以上的验证，即得抽取 h 个成分时第 j 个因变量 y ( j 1,2, , p) j = " 的 预测误差平方和为 ∑= = − n i j h yij y i j h 1 2 ( ) PRESS ( ) ( ˆ ( )) （ j = 1,2,", p ） T p Y ( y , , y ) = 1 " 的预测误差平方和为 ∑= = p i h j h 1 PRESS( ) PRESS ( ) 。 另外，再采用所有的样本点，拟合含h 个成分的回归方程。这时，记第i 个样本点 的预测值为 yˆ (h) ij ，则可以定义 j y 的误差平方和为 ∑= = − n i j h yij yij h 1 2 SS ( ) ( ˆ ( )) 定义Y 的误差平方和为 ∑= = p j h j h 1 SS( ) SS ( ) 当 PRESS(h) 达到最小值时，对应的 h 即为所求的成分个数。通常，总有 PRESS(h) 大于SS(h) ，而SS(h) 则小于SS(h −1)。因此，在提取成分时，总希望比 值 PRESS(h) SS(h −1) 越小越好；一般可设定限制值为 0.05，即当
PRESSO(h)/SS(h-1)≤(1-0.05)2=0952 ,增加成分L有利于模型精度的提高。或者反过来说,当 PRESS(h)/SS(h-1)>0.95 时,就认为增加新的成分b,对减少方程的预测误差无明显的改善作用。 为此,定义交叉有效性为Q2=1- PRESS(h)/SSh-1),这样,在建模的每一步 计算结束前,均进行交叉有效性检验,如果在第h步有Q2<1-0.952=00985,则模 型达到精度要求,可停止提取成分:若Q2≥00975,表示第h步提取的tn成分的边际 贡献显著,应继续第h+1步计算。 §2一种更简洁的计算方法 上节介绍的算法原则和推导过程的思路在目前的文献中是最为常见的。然而,还有 一种更为简洁的计算方法,即直接在E0,…E-1矩阵中提取成分l1…,L1(r≤m)。要 求L能尽可能多地携带X中的信息,同时,t对因变量系统F有最大的解释能力。注 意,无需在F中提取成分得分,这可以使计算过程大为简化,并且对算法结论的解 释也更为方便。 偏最小二乘法的简记算法的步骤如下: (1)求矩阵 EFFE0最大特征值所对应的特征向量1,求得成分1=v1X 计算成分得分向量=E,和残差矩阵E1=E0-ia1,其中a=E/ (2)求矩阵EFFE1最大特征值所对应的特征向量v2,求得成分t2=2X 计算成分得分向量=E,和残差矩阵E2=E1-1a,其中a3=EF/ (r)至第r步,求矩阵EFFE最大特征值所对应的特征向量v,,求得成分 ln=WX,计算成分得分向量,=Ew
-535- 2 2 PRESS(h) SS(h −1) ≤ (1− 0.05) = 0.95 时，增加成分 ht 有利于模型精度的提高。或者反过来说，当 2 PRESS(h) SS(h −1) > 0.95 时，就认为增加新的成分 ht ，对减少方程的预测误差无明显的改善作用。 为此，定义交叉有效性为 1 PRESS( ) SS( 1) 2 Qh = − h h − ，这样，在建模的每一步 计算结束前，均进行交叉有效性检验，如果在第 h 步有 1 0.95 0.0985 2 2 Qh < − = ，则模 型达到精度要求，可停止提取成分；若 0.0975 2 Qh ≥ ，表示第h 步提取的 ht 成分的边际 贡献显著，应继续第 h +1步计算。 §2 一种更简洁的计算方法 上节介绍的算法原则和推导过程的思路在目前的文献中是最为常见的。然而，还有 一种更为简洁的计算方法，即直接在 0 1 , , E " Er− 矩阵中提取成分 r t , ,t 1 " (r ≤ m )。要 求 ht 能尽可能多地携带 X 中的信息，同时， ht 对因变量系统 F0 有最大的解释能力。注 意，无需在 F0 中提取成分得分uh ，这可以使计算过程大为简化，并且对算法结论的解 释也更为方便。 偏最小二乘法的简记算法的步骤如下： （1）求矩阵 E0 F0F0 E0 T T 最大特征值所对应的特征向量 w1，求得成分t w XT 1 = 1 ， 计算成分得分向量 1 0 1 t ˆ = E w ，和残差矩阵 T E E t 1 0 1 1 = − ˆα ，其中 2 1 0 1 1 E t ˆ t ˆ T α = 。 （2）求矩阵 E1 F0F0 E1 T T 最大特征值所对应的特征向量 w2 ，求得成分t w XT 2 = 2 ， 计算成分得分向量 2 1 2 t ˆ = E w ，和残差矩阵 T E E t 2 1 2 2 = − ˆα ，其中 2 2 1 2 2 E t ˆ t ˆ T α = 。 # （r）至第 r 步，求矩阵 −1 0 0 r−1 T T Er F F E 最大特征值所对应的特征向量 wr ，求得成分 t w XT r = r ，计算成分得分向量 r Er wr t 1 ˆ = −
如果根据交叉有效性,确定共抽取r个成分12…,可以得到一个满意的预测模型, 则求F在12…,上的普通最小二乘回归方程为 F=i1B1+…+1B+F 把l=m1x1+…+mxn(k=12,…,r),代入Y=1B1+…+1,B,即得p个 因变量的偏最小二乘回归方程式 y 这里的v满足=Ev,=∏(-w)mn §3案例分析 本节采用兰纳胡德( Finnerud)给出的关于体能训练的数据进行偏最小二乘回归建 模。在这个数据系统中被测的样本点,是某健身俱乐部的20位中年男子。被测变量分 为两组。第一组是身体特征指标X,包括:体重、腰围、脉搏。第二组变量是训练结 果指标Y,包括:单杠、弯曲、跳高。原始数据见表1。 表1体能训练数据 体重(x) 110 38 58 12 105 8 125 11 13 210 13 14 215 202 210 536-
-536- 如果根据交叉有效性，确定共抽取 r 个成分 r t , ,t 1 " 可以得到一个满意的预测模型， 则求 F0 在 r t ˆ , ,t ˆ 1 " 上的普通最小二乘回归方程为 r T r r T F0 = t ˆ 1β1 +"+ t ˆ β + F 把 k k km m t w x w x * 1 * = 1 +"+ （ k = 1,2,",r ），代入 r r Y = t 1β1 +"+ t β ，即得 p 个 因变量的偏最小二乘回归方程式 j j jm m y = a x +"+ a x 1 1 ，（ j = 1,2,",m）。 这里的 * wh 满足 * 0 ˆ h E wh t = ， ∏ − = = − 1 1 * ( ) h j h T wh I wj α j w 。 §3 案例分析 本节采用兰纳胡德（Linnerud）给出的关于体能训练的数据进行偏最小二乘回归建 模。在这个数据系统中被测的样本点，是某健身俱乐部的 20 位中年男子。被测变量分 为两组。第一组是身体特征指标 X ，包括：体重、腰围、脉搏。第二组变量是训练结 果指标Y ，包括：单杠、弯曲、跳高。原始数据见表 1。 表 1 体能训练数据 No 体重(x1) 腰围(x2) 脉搏(x3) 单杠(y1) 弯曲(y2) 跳高(y3) 1 191 36 50 5 162 60 2 189 37 52 2 110 60 3 193 38 58 12 101 101 4 162 35 62 12 105 37 5 189 35 46 13 155 58 6 182 36 56 4 101 42 7 211 38 56 8 101 38 8 167 34 60 6 125 40 9 176 31 74 15 200 40 10 154 33 56 17 251 250 11 169 34 50 17 120 38 12 166 33 52 13 210 115 13 154 34 64 14 215 105 14 247 46 50 1 50 50 15 193 36 46 6 70 31 16 202 37 62 12 210 120
17 176 37 54 60 52 均值 35.4 145.55 标准差 4.6905 3.202 72104 5.2863 62.5666 51.2775 表2给出了这6个变量的简单相关系数矩阵。从相关系数矩阵可以看出,体重与腰 围是正相关的:体重、腰围与脉搏负相关;而在单杠、弯曲与跳高之间是正相关的。从 两组变量间的关系看,单杠、弯曲和跳高的训练成绩与体重、腰围负相关,与脉搏正相 表2相关系数矩阵 702-0.3658-0.3897-0.4931-0.2263 0.3529-0.5522-0.6456-0.1915 0.3658 0.3529 0.1506 0.225 0.0349 0.3897 0.1506 0.6957 0.4958 0.4931 0.6456 0.6957 0.6692 0.2263 0.1915 0.0349 0.4958 0.6692 利用如下的 MATLAB程序: load pz. txt%原始数据存放在纯文本文件 pz. txt中 mu=mean(pz),sig=std(pz),%求均值和标准差 rr=corrcoef(pz),%求相关系数矩阵 data= Ascore(pz),%数据标准化 n=3,m=3;,%n是自变量的个数,m是因变量的个数 x0=pz(, 1: n): yo=pz( n +1: end); um=size(e0,1)%求样本点的个数 chg=eye(n),%w到w*变换矩阵的初始化 for F=l: n %以下计算w,w和t的得分向量 matrix=e0*f0*fo*e0 Ivec, vall==eig( matrIX,%求特征值和特征向量 val=diag(va),%提出对角线元素 Ival, ind] =sort(val, ' descend); w(i)=vecc;ind(1),%提出最大特征值对应的特征向量
-537- 17 176 37 54 4 60 25 18 157 32 52 11 230 80 19 156 33 54 15 225 73 20 138 33 68 2 110 43 均值 标准差 178.6 24.6905 35.4 3.202 56.1 7.2104 9.45 5.2863 145.55 62.5666 70.3 51.2775 表 2 给出了这 6 个变量的简单相关系数矩阵。从相关系数矩阵可以看出，体重与腰 围是正相关的；体重、腰围与脉搏负相关；而在单杠、弯曲与跳高之间是正相关的。从 两组变量间的关系看，单杠、弯曲和跳高的训练成绩与体重、腰围负相关，与脉搏正相 关。 表 2 相关系数矩阵 1 0.8702 -0.3658 -0.3897 -0.4931 -0.2263 0.8702 1 -0.3529 -0.5522 -0.6456 -0.1915 -0.3658 -0.3529 1 0.1506 0.225 0.0349 -0.3897 -0.5522 0.1506 1 0.6957 0.4958 -0.4931 -0.6456 0.225 0.6957 1 0.6692 -0.2263 -0.1915 0.0349 0.4958 0.6692 1 利用如下的 MATLAB 程序： clc,clear load pz.txt %原始数据存放在纯文本文件 pz.txt 中 mu=mean(pz);sig=std(pz); %求均值和标准差 rr=corrcoef(pz); %求相关系数矩阵 data=zscore(pz); %数据标准化 n=3;m=3; %n 是自变量的个数,m 是因变量的个数 x0=pz(:,1:n);y0=pz(:,n+1:end); e0=data(:,1:n);f0=data(:,n+1:end); num=size(e0,1);%求样本点的个数 chg=eye(n); %w 到 w*变换矩阵的初始化 for i=1:n %以下计算 w，w*和 t 的得分向量， matrix=e0'*f0*f0'*e0; [vec,val]=eig(matrix); %求特征值和特征向量 val=diag(val); %提出对角线元素 [val,ind]=sort(val,'descend'); w(:,i)=vec(:,ind(1)); %提出最大特征值对应的特征向量
w star(:,i)=chg*w(i),%计算w*的取值 t(:,i=e0*w(:,i) %计算成分ti的得分 alpha=e0*(;)/(,y*t(:,1);,%计算 alpha_i chg=chg (eye(n)-w( 1)*alpha %计算w到w*的变换矩阵 e=e0-t(:)* alpha,%计算残差矩阵 %以下计算sSi的值 beta=[,1)ones(mum,1)f0,%求回归方程的系数 beta(end,-}=:%删除回归分析的常数项 cancha=f0-t(,l;i)*beta;%求残差矩阵 sS=sum( sum(cancha^2);%求误差平方和 %以下计算 press( for j=l: num tl=t(:,1:1);fl=f0 shet=t(j;),shef=f(G,),%把舍去的第j个样本点保存起来 t()=f(-)= %删除第j个观测值 betal- Itl, ones(num-1,1);%求回归分析的系数 betal(end, 0: %删除回归分析的常数项 cancha= she f-she t*betal;%求残差向量 press IOF press(isum(press 1) Q h2(iFl-press(i)ss(i-1); Qh2(1)=1 fQh2()<0.0975 fprintf提出的成分个数r=%d,i); beta z=t(:,1r)ones(mum,1)t0,%求Y关于t的回归系数 beta_z(end,)=[%删除常数项 xishu= w star(,l;r)* beta z,%求Y关于X的回归系数,且是针对标准数据的回归系数 每一列是一个回归方程 mu x=mu(l: n) mu y=mu(n+1: end) sig x=sig(l: n), sig y=sig(n+l: end) 538
-538- w_star(:,i)=chg*w(:,i); %计算 w*的取值 t(:,i)=e0*w(:,i); %计算成分 ti 的得分 alpha=e0'*t(:,i)/(t(:,i)'*t(:,i)); %计算 alpha_i chg=chg*(eye(n)-w(:,i)*alpha'); %计算 w 到 w*的变换矩阵 e=e0-t(:,i)*alpha'; %计算残差矩阵 e0=e; %以下计算 ss(i)的值 beta=[t(:,1:i),ones(num,1)]\f0; %求回归方程的系数 beta(end,:)=[]; %删除回归分析的常数项 cancha=f0-t(:,1:i)*beta; %求残差矩阵 ss(i)=sum(sum(cancha.^2)); %求误差平方和 %以下计算 press(i) for j=1:num t1=t(:,1:i);f1=f0; she_t=t1(j,:);she_f=f1(j,:); %把舍去的第 j 个样本点保存起来 t1(j,:)=[];f1(j,:)=[]; %删除第 j 个观测值 beta1=[t1,ones(num-1,1)]\f1; %求回归分析的系数 beta1(end,:)=[]; %删除回归分析的常数项 cancha=she_f-she_t*beta1; %求残差向量 press_i(j)=sum(cancha.^2); end press(i)=sum(press_i); if i>1 Q_h2(i)=1-press(i)/ss(i-1); else Q_h2(1)=1; end if Q_h2(i)<0.0975 fprintf('提出的成分个数 r=%d',i); r=i; break end end beta_z=[t(:,1:r),ones(num,1)]\f0; %求 Y 关于 t 的回归系数 beta_z(end,:)=[]; %删除常数项 xishu=w_star(:,1:r)*beta_z; %求Y关于X的回归系数，且是针对标准数据的回归系数， 每一列是一个回归方程 mu_x=mu(1:n);mu_y=mu(n+1:end); sig_x=sig(1:n);sig_y=sig(n+1:end);
for F =l ch(=muy(i)-mux/igx*sigy()* xishu(n),%计算原始数据的回归方程的常数 项 xish(:,i= xishu(,i)/sgx*sigy(i),%计算原始数据的回归方程的系数,每一列是一 个回归方程 l=ch0xish]%显示回归方程的系数,每一列是一个方程,每一列的第一个数是 常数项 save mydata xo yO num xishu ch0 xish 计算得只要提出两个成分,l2即可,交叉有效性Q2=-0.969。Wn与v2的取值 见表3,成分t的得分b见表4 表3h与h的取值 自变量 0.5899 -04688 -0.5899 -0.3679 0.7713 0.5680 0.7713 0.6999 0.6765 0.2389 0.6356 表4成分L的得分t -0.6429-0.7697-0.90740.6884-0.4867-.2211-1.40370.74361.71511.1626 2-0.5914-0.16670.52120.68-.13280.07170.07670.21060.6549-0.1668 11 12 13 14 16 17 18 19 20 i1|0.36450.74331.1867-4398-0.8232-0.749-0.39291.1991.04851.924 -0.7007-0.6830.75 0.76 -0.97380.52110.2034-0.7827-0.37291.1294 标准化变量关于成分t1的回归模型如下 D=n 4,+r2k42, k
-539- for i=1:m ch0(i)=mu_y(i)-mu_x./sig_x*sig_y(i)*xishu(:,i); %计算原始数据的回归方程的常数 项 end for i=1:m xish(:,i)=xishu(:,i)./sig_x'*sig_y(i); %计算原始数据的回归方程的系数，每一列是一 个回归方程 end sol=[ch0;xish] %显示回归方程的系数，每一列是一个方程，每一列的第一个数是 常数项 save mydata x0 y0 num xishu ch0 xish 计算得只要提出两个成分 1 2 t ,t 即可，交叉有效性 0.1969 2 Q2 = − 。 wh 与 * wh 的取值 见表 3，成分 ht 的得分 ht ˆ 见表 4。 表 3 wh 与 * wh 的取值 自变量 w1 w2 * 1 w * 2 w x1 -0.5899 -0.4688 -0.5899 -0.3679 x2 -0.7713 0.5680 -0.7713 0.6999 x3 0.2389 0.6765 0.2389 0.6356 表 4 成分 ht 的得分 ht ˆ No 1 2 3 4 5 6 7 8 9 10 1 t ˆ -0.6429 -0.7697 -0.9074 0.6884 -0.4867 -0.2291 -1.4037 0.7436 1.7151 1.1626 2 t ˆ -0.5914 -0.1667 0.5212 0.68 -1.1328 0.0717 0.0767 0.2106 0.6549 -0.1668 No 11 12 13 14 15 16 17 18 19 20 1 t ˆ 0.3645 0.7433 1.1867 -4.3898 -0.8232 -0.749 -0.3929 1.1993 1.0485 1.9424 2 t ˆ -0.7007 -0.6983 0.757 0.76 -0.9738 0.5211 0.2034 -0.7827 -0.3729 1.1294 标准化变量 k y ~ 关于成分 1t 的回归模型如下 1 1 2 2 ~y r t r t k = k + k ， k = 1,2,3
由于成分t可以写成原变量的标准化变量x的函数,即有 Ih=wwx, t w2hx2+ W3hx3 由此可得由成分1所建立的偏最小二乘回归模型为 丙=(11+21x2+31x3)+24(2x+22+32x3) =(r11+n242)x1+(r121+E212)x2+(rn2+E4n2)3 有关h1=(rh22,3)的计算结果见表5。 表5回归系数r 0.3416 0.4161 0.1430 0.3364 -0.2908 0.0652 所以,有 y1=-0.0778x1-0.4989x2-0.1322x3 y2=-0.1385x1-0.524x2-0.0854x3 3=-0.0604x1-0.1559x2-0.0073x3 将标准化变量y,x(k=12,3)分别还原成原始变量y,x4(k=1,2,3),则回归方 程为 y1=470197-0.0167x1-0.8237x2-0.0969x3 y2=6125671-0.3509x1-10.2477x2-0.7412x3 y3=1839849-0.1253x-2.4969x2-00518x 为了更直观、迅速地观察各个自变量在解释y4(k=12,3)时的边际作用,可以绘 制回归系数图,见图1。这个图是针对标准化数据的回归方程的。 从回归系数图中可以立刻观察到,腰围变量在解释三个回归方程时起到了极为重要 的作用。然而,与单杠及弯曲相比,跳高成绩的回归方程显然不够理想,三个自变量对 它的解释能力均很低。 5
-540- 由于成分 ht 可以写成原变量的标准化变量 j x ~ 的函数，即有 3 * 2 3 * 1 2 * 1 ~ ~ ~ t w x w x w x h = h + h + h 由此可得由成分 1t 所建立的偏最小二乘回归模型为 3 * 2 32 * 2 1 31 * 2 22 * 1 1 21 * 2 12 * 1 11 3 * 2 32 * 1 22 * 3 2 12 * 2 31 * 1 21 * 1 11 ~ ( ) ~ ( ) ~ ( ) ) ~ ~ ~ ) ( ~ ~ ~ ( ~ r w r w x r w r w x r w r w x y r w x w x w x r w x w x w x k k k k k k k k k = + + + + + = + + + + + 有关 ( , , ) h h1 h2 h3 r = r r r 的计算结果见表 5。 表 5 回归系数 hr k 1 2 3 r1 0.3416 0.4161 0.1430 r2 -0.3364 -0.2908 -0.0652 所以，有 1 1 2 3 ~ 0.1322 ~ 0.4989 ~ 0.0778 ~y = − x − x − x 2 1 2 3 ~ 0.0854 ~ 0.5244 ~ 0.1385 ~y = − x − x − x 3 1 2 3 ~ 0.0073 ~ 0.1559 ~ 0.0604 ~y = − x − x − x 将标准化变量 ( 1,2,3) ~, ~yk xk k = 分别还原成原始变量 y , x (k = 1,2,3) k k ，则回归方 程为 1 1 2 0969 3 y = 47.0197 − 0.0167x − 0.8237x − 0. x 2 1 2 7412 3 y = 612.5671− 0.3509x −10.2477x − 0. x 3 1 2 0518 3 y =183.9849 − 0.1253x − 2.4969x − 0. x 为了更直观、迅速地观察各个自变量在解释 y (k =1,2,3) k 时的边际作用，可以绘 制回归系数图，见图 1。这个图是针对标准化数据的回归方程的。 从回归系数图中可以立刻观察到，腰围变量在解释三个回归方程时起到了极为重要 的作用。然而，与单杠及弯曲相比，跳高成绩的回归方程显然不够理想，三个自变量对 它的解释能力均很低
