原理

偏最小二乘回歸的思想是將 Y = [ y₁, y₂, ... , y_n ] 椎椰，X = [ x₁, x₂, ... , x_n ]

步驟

1. 將 X 和 Y 進(jìn)行標(biāo)準(zhǔn)化

設(shè) X 標(biāo)準(zhǔn)化后的矩陣為 E₀痹换， Y 標(biāo)準(zhǔn)化后的矩陣為 F₀

2. 求解第一個(gè)成分

設(shè) E₀ 和 F₀ 的第一成分為 t₁ 和 u₁。這一步的目的是使得 E₀ 和 F₀ 的第一成分之間的相關(guān)性達(dá)到最大嗤练，方便用 E₀ 的第一成分 t₁ 線性表示 F₀

那么 E₀ 和 F₀ 的第一個(gè)成分 t₁胖眷，u₁ 定義如下：

3. 建立第一個(gè)成分與 E₀ 和 F₀ 之間的關(guān)系

根據(jù)最小二乘法原理求解系數(shù)得：

4. 計(jì)算殘差矩陣

由式子 5-11 有：

一般地有

其中：

1. E₀ 代表決策變量 X 標(biāo)準(zhǔn)化后的矩陣竭恬，F(xiàn)₀ 代表響應(yīng)變量 F 標(biāo)準(zhǔn)化后的矩陣
2. E_1-i袱箱，F(xiàn)_1-i 代表殘差矩陣境析，參考式子 5-13 到 5-16

5. 重構(gòu) X，Y 和成分 t 的關(guān)系

以下的 X咆槽，Y 原始數(shù)據(jù)矩陣默認(rèn)為標(biāo)準(zhǔn)化后的：
因此經(jīng)過多次分解以后：

6. 重構(gòu) X 與 Y 之間的關(guān)系

以下的 X陈轿，Y 原始數(shù)據(jù)矩陣默認(rèn)為標(biāo)準(zhǔn)化后的：
首先理解偏最小二乘的外部結(jié)構(gòu)和內(nèi)部結(jié)構(gòu)：

假設(shè)每個(gè)成分 t 與 X 的關(guān)系為：

即：

整理得到 X 與 Y 之間的關(guān)系

得證

理解 t 與 u 之間的關(guān)系

t 與 u 之間滿足相關(guān)性最大，這樣可以用 t 替代 u 來線性表示 Y

代碼

import numpy as np

# 設(shè)決策變量 X = { x1, x2, x3 }罗晕，響應(yīng)變量 Y = { y1, y2, y3 }
x1=[191,189,193,162,189,182,211,167,176,154,169,166,154,247,193,202,176,157,156,138]
x2=[36,37,38,35,35,36,38,34,31,33,34,33,34,46,36,37,37,32,33,33]
x3=[50,52,58,62,46,56,56,60,74,56,50,52,64,50,46,62,54,52,54,68]
y1=[5,2,12,12,13,4,8,6,15,17,17,13,14,1,6,12,4,11,15,2]
y2=[162,110,101,105,155,101,101,125,200,251,120,210,215,50,70,210,60,230,225,110]
y3=[60,60,101,37,58,42,38,40,40,250,38,115,105,50,31,120,25,80,73,43]


data_raw=np.array([x1,x2,x3,y1,y2,y3])
data_raw=data_raw.T # 輸入原始數(shù)據(jù)济欢，行數(shù)為樣本數(shù)，列數(shù)為特征數(shù)
# 數(shù)據(jù)標(biāo)準(zhǔn)化
num=np.size(data_raw,0) # 樣本個(gè)數(shù)
mu=np.mean(data_raw,axis=0) # 按列求均值
sig=(np.std(data_raw,axis=0)) # 按列求標(biāo)準(zhǔn)差
data=(data_raw-mu)/sig # 標(biāo)準(zhǔn)化小渊，按列減去均值除以標(biāo)準(zhǔn)差
# 提取自變量和因變量數(shù)據(jù)

n=3 # 自變量個(gè)數(shù)
m=3 # 因變量個(gè)數(shù)
x0=data_raw[:,0:n] # 原始的自變量數(shù)據(jù)
y0=data_raw[:,n:n+m] # 原始的變量數(shù)據(jù)
e0=data[:,0:n] # 標(biāo)準(zhǔn)化后的自變量數(shù)據(jù)
f0=data[:,n:n+m] # 標(biāo)準(zhǔn)化后的因變量數(shù)據(jù)

chg=np.eye(n) # w 到 w* 變換矩陣的初始化

w=np.empty((n,0)) # 初始化投影軸矩陣
w_star=np.empty((n, 0)) # w* 矩陣初始化
t=np.empty((num, 0)) # 得分矩陣初始化
ss=np.empty(0) # 或者ss=[]，誤差平方和初始化
press=[] # 預(yù)測(cè)誤差平方和初始化
Q_h2=np.zeros(n) # 有效性判斷條件值初始化
print(Q_h2)

for i in range(n):
    matrix=e0.T@f0@f0.T@e0 # 構(gòu)造矩陣E'FF'E
    val,vec=np.linalg.eig(matrix) # 計(jì)算特征值和特征向量
    index=np.argsort(val)[::-1] # 獲取特征值從大到小排序前的索引
    val_sort=val[index] # 特征值由大到小排序
    vec_sort=vec[:,index] # 特征向量按照特征值的順序排列
    w=np.append(w,vec_sort[:,0][:,np.newaxis],axis=1) # 儲(chǔ)存最大特征向量
    w_star=np.append(w_star,chg@w[:,i][:,np.newaxis],axis=1) # 計(jì)算 w* 的取值
    t=np.append(t,e0@w[:,i][:,np.newaxis],axis=1) # 計(jì)算投影
    alpha=e0.T@t[:,i][:,np.newaxis]/(t[:,i]@t[:,i]) # 計(jì)算自變量和主成分之間的回歸系數(shù)
    chg=chg@(np.eye(n)-(w[:,i][:,np.newaxis]@alpha.T)) # 計(jì)算 w 到 w*的變換矩陣
    e1=e0-t[:,i][:,np.newaxis]@alpha.T # 計(jì)算殘差矩陣
    e0=e1 # 更新殘差矩陣


    beta=np.linalg.pinv(t)@f0 #求回歸方程的系數(shù)茫叭，數(shù)據(jù)標(biāo)準(zhǔn)化酬屉，沒有常數(shù)項(xiàng)
    res=np.array(f0-t@beta) #求殘差
    ss=np.append(ss,np.sum(res**2))#殘差平方和
    #-----求解殘差平方和press
    press_i=[] #初始化誤差平方和矩陣
    for j in range(num):
        t_inter=t[:,0:i+1]
        f_inter=f0
        t_inter_del=t_inter[j,:] #把舍去的第 j 個(gè)樣本點(diǎn)保存起來,自變量
        f_inter_del=f_inter[j,:] #把舍去的第 j 個(gè)樣本點(diǎn)保存起來，因變量
        t_inter= np.delete(t_inter,j,axis=0) #刪除自變量第 j 個(gè)觀測(cè)值
        f_inter= np.delete(f_inter,j,axis=0) #刪除因變量第 j 個(gè)觀測(cè)值
        t_inter=np.append(t_inter,np.ones((num-1,1)),axis=1)
        beta1=np.linalg.pinv(t_inter)@f_inter # 求回歸分析的系數(shù),這里帶有常數(shù)項(xiàng)
        res=f_inter_del-t_inter_del[:,np.newaxis].T@beta1[0:len(beta1)-1,:]-beta1[len(beta1)-1,:] #計(jì)算殘差
        res=np.array(res)
        press_i.append(np.sum(res**2)) #殘差平方和揍愁，并存儲(chǔ)
    press.append(np.sum(press_i)) #預(yù)測(cè)誤差平方和

    Q_h2[0]=1
    if i>0:
        Q_h2[i]=1-press[i]/ss[i-1]
    if Q_h2[i]<0.0975:
        print('提出的成分個(gè)數(shù) r=',i+1)
        break


beta_Y_t=np.linalg.pinv(t)@f0 #求Y*關(guān)于t的回歸系數(shù)
beta_Y_X=w_star@beta_Y_t#求Y*關(guān)于X*的回歸系數(shù)
mu_x=mu[0:n] #提取自變量的均值
mu_y=mu[n:n+m] #提取因變量的均值
sig_x=sig[0:n] #提取自變量的標(biāo)準(zhǔn)差
sig_y=sig[n:n+m] #提取因變量的標(biāo)準(zhǔn)差
ch0=mu_y-mu_x[:,np.newaxis].T/sig_x[:,np.newaxis].T@beta_Y_X*sig_y[:,np.newaxis].T#算原始數(shù)據(jù)回歸方程的常數(shù)項(xiàng)
beta_target=np.empty((n,0)) #回歸方程的系數(shù)矩陣初始化

for i in range(m):
    a=beta_Y_X[:,i][:,np.newaxis]/sig_x[:,np.newaxis]*sig_y[i]#計(jì)算原始數(shù)據(jù)回歸方程的系數(shù)
    beta_target=np.append(beta_target,a,axis=1)
target=np.concatenate([ch0,beta_target],axis=0) #回歸方程的系數(shù)呐萨，每一列是一個(gè)方程，每一列的第一個(gè)數(shù)是常數(shù)項(xiàng)

print(target)