在线时间:8:00-16:00
迪恩网络APP
随时随地掌握行业动态
扫描二维码
关注迪恩网络微信公众号
本文所有测试用数据集均来自Matlab,并使用Matlab封装的回归方法,对自己实现的代码做了验证,本文参考文献及资料如下: Reference: [1] GELADI P, KOWALSKI B R. Partial least-squares regression: a tutorial [J]. Analytica chimica acta, 1986, 185(1-17). [2] WU F Y, ASADA H H. Implicit and intuitive grasp posture control for wearable robotic fingers: A data-driven method using partial least squares [J]. IEEE Transactions on Robotics, 2016, 32(1): 176-86. [3] https://en.wikipedia.org/wiki/Ordinary_least_squares [4] https://en.wikipedia.org/wiki/Principal_component_analysis [5] https://en.wikipedia.org/wiki/Principal_component_regression [6] https://en.wikipedia.org/wiki/Partial_least_squares_regression
完整Matlab代码实现: https://github.com/ShieldQiQi/PCA-PCR-PLSR-Matlab-code 一、OLSR 即为普通最小二乘回归,对此我们应该十分熟悉,各种大物材料力学实验都会用到这种方法,只不过我们当时使用的单变量的数据,当数据集涉及到矩阵,多维变量的形式时,就需要使用更加普遍适用的模型,我们设原始数据自变量(independent value)矩阵为$ X∈R_{n{\times}m} $,即X数据集含有n个样本,每个样本有m个特征变量;设原始数据因变量(dependent value)矩阵为$ Y∈R_{n{\times}p} $,即Y数据集含有n个样本,每个样本有p个特征变量。构建的最小二乘回归模型为: $$ Y=X{\cdot}B+E \tag{1} $$ 上式中$ B∈R_{m{\times}p} $为回归模型的系数矩阵,$ E∈R_{n{\times}p} $为模型预测的残差。B的通用解法参考维基百科,为: $$ B=(X^{T}X)^{-1}X^{T}Y \tag{2} $$ 二、PCA PCA本质上是一种建立一种维度小于原始数据维度(特征变量数)正交基底空间,将原始数据投影到新的低维空间,以达到数据降维而保持原有特性的方法。PCA的步骤为:
至于为什么这样做,PCA的原理可以参考维基百科,或者我的这篇博文: https://www.cnblogs.com/QiQi-Robotics/p/14303718.html 在实际应用中,计算协方差矩阵的特征向量常采用迭代计算的方式,常用的方法为NIPALS,Matlab精简代码(Matlab使用的为散布矩阵,而我的代码为协方差矩阵,所以特征值会相差(n-1)倍)实现如下: 1 % 迭代得到num个成份 2 for h = 1:num 3 % step(1) 4 % --------------------------------------------------------------------- 5 % 取T(:,h)为任意一个X_centered中的列向量,此处直接取第一列 6 T(:,h) = X_iteration(:,1); 7 8 % step(2) to step(5) 9 % 迭代直到收敛到容忍度内的主成分 10 while(1) 11 P(:,h) = X_iteration'*T(:,h)/(T(:,h)'*T(:,h)); 12 % 归一化P(:,h) 13 P(:,h) = P(:,h)/sqrt(P(:,h)'*P(:,h)); 14 t_temp = T(:,h); 15 T(:,h) = X_iteration*P(:,h)/(P(:,h)'*P(:,h)); 16 17 % 检查当前T(:,h)与上一步T(:,h)是否相等以决定是否继续迭代 18 if max(abs(T(:,h)-t_temp)) <= tolerance 19 % 存储按顺序排列的特征值 20 % 注意此处的特征值为协方差矩阵的特征值,而matlab PCA方法使用的为散布矩阵(离散度矩阵),故后者的特征值为前者的(n-1)倍 21 eigenValues(h) = P(:,h)'*(X_centered'*X_centered)*P(:,h); 22 break; 23 else 24 end 25 end 26 27 % 计算残差,更新数据矩阵 28 % --------------------------------------------------------------------- 29 X_iteration = X_iteration - T(:,h)*P(:,h)'; 30 end 三、PCR PCR使用的回归方法是OLSR,只不过回归的模型是建立在主成分空间,以防止原始数据的共线性问题导致模型建立不准确,步骤如下:
Matlab精简代码如下: 1 % 定义测试集样本的数量 2 r = n; 3 % 将原始数据降维到主成分空间(T)后,使用OLS最小二乘回归获取系数矩阵 4 B_inPca = inv(T'*T)*T'*Y_centered; 5 %B_inPca = regress(Y-mean(Y), T(:,1:num)); 6 % 将系数矩阵从主成分空间转化到原始空间 7 B_estimated = P*B_inPca; 8 9 % 定义测试集,此处直接使用原始数据的前r行 10 X_validate = zeros(r,m); 11 % 对原始数据集居中列平均化 12 for j = 1:m 13 % 注意,此处减去的平均值应该为模型数据集的平均值,而非新数据的平均值 14 X_validate(:,j) = X(1:r,j) - mean(X(:,j)); 15 end 16 17 Y_estimated = X_validate*B_estimated; 18 for i = 1:p 19 % 注意此处最终的输出需要加上数据集Y的均值 20 Y_estimated(:,i) = Y_estimated(:,i) + mean(Y(:,i)); 21 end 四、PLSR PLSR相对于PCR的一个优点在于在使用更少的主成分可以获得更具有鲁棒性的预测结果(具体可以查看Matlab中关于PLSR的帮助文档),具体步骤查阅论文 [1]。精简版Matlab代码如下: 1.建立模型部分 1 % 迭代得到num个成份 2 for h = 1:num 3 % step(1) 4 % --------------------------------------------------------------------- 5 % 取u_h为任意一个Y_centered中的列向量,此处直接取第一列 6 U(:,h) = Y_centered(:,1); 7 8 % step(2) to step(8) 9 % --------------------------------------------------------------------- 10 while 1 11 % 在数据矩阵X_centered中 12 W(:,h) = X_centered'*U(:,h)/(U(:,h)'*U(:,h)); 13 % 对数据进行归一化 14 W(:,h) = W(:,h)/sqrt(W(:,h)'*W(:,h)); 15 t_temp = T(:,h); 16 T(:,h) = X_centered*W(:,h)/(W(:,h)'*W(:,h)); 17 18 % 在数据矩阵Y_centered中 19 Q(:,h) = Y_centered'*T(:,h)/(T(:,h)'*T(:,h)); 20 % 对数据进行归一化 21 Q(:,h) = Q(:,h)/sqrt(Q(:,h)'*Q(:,h)); 22 U(:,h) = Y_centered*Q(:,h)/(Q(:,h)'*Q(:,h)); 23 24 % 检查T(:,h)与T(:,h)的前一步是否相等,若小于某个数值则该PLS成份迭代完成,否则返回继续迭代 25 if max(abs(T(:,h)-t_temp)) <= tolerance 26 break; 27 else 28 end 29 end 30 31 % step(9) to step(13) 32 % --------------------------------------------------------------------- 33 P(:,h) = X_centered'*T(:,h)/(T(:,h)'*T(:,h)); 34 % 对数据进行归一化 35 p_norm = sqrt(P(:,h)'*P(:,h)); 36 P(:,h) = P(:,h)/p_norm; 37 T(:,h) = T(:,h)*p_norm; 38 W(:,h) = W(:,h)*p_norm; 39 B(h) = U(:,h)'*T(:,h)/(T(:,h)'*T(:,h)); 40 41 % 计算残差,更新数据矩阵 42 % --------------------------------------------------------------------- 43 X_centered = X_centered - T(:,h)*P(:,h)'; 44 Y_centered = Y_centered - B(h)*T(:,h)*Q(:,h)'; 45 end 2.预测部分 1 % 对原始数据集居中列平均化 2 for j = 1:m 3 % 注意,此处减去的平均值应该为模型数据集的平均值,而非新数据的平均值 4 X_validate(1:r,j) = X(1:r,j) - mean(X(:,j)); 5 end 6 7 % 计算预测的T 8 for h = 1:num 9 T_est(:,h) = X_validate*W(:,h); 10 X_validate = X_validate - T_est(:,h)*P(:,h)'; 11 end 12 13 % 计算预测的Y 14 for h = 1:num 15 Y_estimated = Y_estimated + B(h)*T_est(:,h)*Q(:,h)'; 16 end 17 for i = 1:p 18 % 注意此处最终的输出需要加上数据集Y的均值 19 Y_estimated(:,i) = Y_estimated(:,i) + mean(Y(:,i)); 20 end 五、实验结果
图1 Matlab PLSR算法(SIMPLS)和自定义PLSR(NIPALS)方法效果对比
|
2023-10-27
2022-08-15
2022-08-17
2022-09-23
2022-08-13
请发表评论