【Matlab学习手记】偏最小二乘回归

偏最小二乘程序，可移植性强。

function [T, P, U, Q, B, W] = PLS(X, Y)
% 函数功能：偏最小二乘
% ===============================================================
% 输入：
%   X：自变量；
%   Y：因变量；
%   TolY：收敛条件，默认为1e-10；
% 输出：
%   T：X的得分矩阵；
%   P：X的负荷矩阵；
%   U：Y的得分矩阵；
%   Q：Y的负荷矩阵；
%   B：回归系数矩阵；
%   W：X的权重矩阵；
% 根据PLS模型预测：Y = X*(P*B*Q')
% 调用示例：
%{
clear;clc;
x = [4 9 6 7 7 8 3 2;6 15 10 15 17 22 9 4;8 21 14 23 27 36 15 6;
10 21 14 13 11 10 3 4; 12 27 18 21 21 24 9 6; 14 33 22 29 31 38 15 8;
16 33 22 19 15 12 3 6; 18 39 26 27 25 26 9 8;20 45 30 35 35 40 15 10];
y = [1 1;3 1;5 1;1 3;3 3;5 3;1 5;3 5;5 5];
N = size(x,1);
% 标准化
xmean = mean(x);
xstd = std(x);
ymean = mean(y);
ystd = std(y);
X = (x - repmat(xmean, N, 1)) ./ repmat(xstd, N, 1);
Y = (y - repmat(ymean, N, 1)) ./ repmat(ystd, N, 1);
[T, P, U, Q, B, W] = PLS(X, Y);
% 拟合值
yp = (x - repmat(xmean, N, 1)) ./ repmat(xstd, N, 1) * (P*B*Q');
% 还原
yp = yp.*repmat(ystd, N, 1) + repmat(ymean, N, 1);
bar([y(:, 1), yp(:, 1)]);
%}
% ===============================================================
TolY = 1e-10;
[rX, cX] = size(X);
[rY, cY] = size(Y);
assert(rX == rY, 'X的行数与Y的行数不一致！');
%%%
n = max(cX, cY);
T = zeros(rX, n);
P = zeros(cX, n);
U = zeros(rY, n);
Q = zeros(cY, n);
B = zeros(n, n);
W = P;
k = 0;
TolX = 1e-10;
% iteration loop if residual is larger than specfied
while norm(Y) > TolY && k < n
    % choose the column of x has the largest square of sum as t.
    % choose the column of y has the largest square of sum as u.    
    [~, tidx] =  max(sum(X.*X));
    [~, uidx] =  max(sum(Y.*Y));
    t1 = X(:, tidx);
    u = Y(:, uidx);
    t = t1 - t1;
    % iteration for outer modeling until convergence
    while norm(t1 - t) > TolX
        w = X'*u;
        w = w/norm(w);
        t = t1;
        t1 = X*w;
        q = Y'*t1;
        q = q/norm(q);
        u = Y*q;
    end
    % update p based on t
    t = t1;
    p = X'*t/(t'*t);
    pnorm = norm(p);
    p = p/pnorm;
    t = t*pnorm;
    w = w*pnorm;
    % regression and residuals
    b = u'*t/(t'*t);
    X = X - t*p';
    Y = Y - b*t*q';
    % save iteration results to outputs:
    k = k + 1;
    T(:, k) = t;
    P(:, k) = p;
    U(:, k) = u;
    Q(:, k) = q;
    W(:, k) = w;
    B(k, k) = b;
end
T(:, k + 1:end) = [];
P(:, k + 1:end) = [];
U(:, k + 1:end) = [];
Q(:, k + 1:end) = [];
W(:, k + 1:end) = [];
B = B(1:k, 1:k);