pls1.m

pls for matlab。matlab中实现偏最小二乘法

% 单因变量偏最小二乘回归分析
% 输入:
% x: 自变量矩阵,
% y: 因变量列向量,如果y是矩阵, 则取第一列作因变量列向量
% step: 提取的主成分数
% 输出:
% A: 回归方程系数, 常数项在前
% a: 标准化数据的回归方程系数
% T: 主成分

% 
% simple partial least squares regression, 
% INPUT:
% x: matrix of independent variables, 
% y: column vector of dependent variable, if y is matrix, take the first column as dependent variable.
% step: component number
% OUTPUT:
% A:  regression coefficent vector for original data
% a:  regression coefficent vector for standardized data
% T:  matrix of principle components

% 
% 
% 
function [A, a, T] = pls1 (x,y,step)

%---------------------------------------------------------------

if nargin ~= 2 && nargin ~= 3
    disp('错误: 参数不正确!')
    return
end

[n,m]= size(x);
[yn,ym] = size(y);

if nargin == 3 && ( step > m || step < 1)
     disp('错误: 参数不正确!')
     return;
elseif nargin == 2
    step = m;
end

if n ~= yn
    disp('错误: 自变量与因变量样本数不一致!')
    return
elseif m<2
    disp('错误: 自变量个数太少,不能少于2!')
    return
elseif n<2
    disp('错误: 样本数太少,不能少于2!')
    return
end
%-------------------------------------------------------------
aver = mean(x);
stdcov = std(x); 
E0 = ( x - aver(ones(n,1),:) )./ stdcov( ones(n,1),:);   % 标准化自变量

if ym == 1
    y1 = y;
else
    y1 = y(:,1);
end

F0 = (y1 - mean(y1) )./std(y1);  % 标准化因变量

%--------------------------------------------------------------
E = E0;
F = F0;

for h = 1:step
    w = ( E'* F )/ norm( E'*F );
    W(:, h ) = w;
    
    t = E * w;
    T(:, h ) = t;     % 主成分
    
    p = E'* t / ( norm( t)^2 );
    P(:, h ) = p;
    
    r = F'* t / ( norm( t)^2 );
    R( h ) = r;       % 系数
    
    E = E - t * p';
    F = F - t * r;
    
    if step == m && nargin == 2 
        if h == step       % 已经是第step步了, 不用再做显著性检验了
            break;
        end
        
        SSi = sum( (F0 - T*R').^2 );
        PRESSi1 = GetPRESS(E0,F0,h+1);
        
        %   PRESSi1                               % 打出来看看
           Q2 = ( 1- PRESSi1/SSi)                % 打出来看看
        
        if ( 1- PRESSi1/SSi) < 0.0975  % 交叉显著性检验
            break;                  % 不显著
        end       
    end
end

% T             % 打出来看看

% 求W* 的值
for i = 1: h   
    W_h_s = eye(m);
    for j = 1:(i-1)
        W_h_s = W_h_s * ( eye(m) - W(:,j) *( P(:,j) )');
    end
    
    W_h_s = W_h_s * W(:,i);  % Wh* 的值
    W_s(:,i) = W_h_s;
end

a_s =  R * W_s'; % 标准化的回归方程系数

a = a_s;              % 返回值

ai = a_s * std(y1) ./ stdcov;
a0 = mean(y1) - sum( ( a_s * std(y1) ./ stdcov ) .* mean(x) ); % 常数

A = [a0,ai];    % 返回值


function re = GetPRESS(E0,F0,h)

[n,m] = size(E0);
F0n = size(F0);
if  F0n ~= n | n <2 | m < 2
    disp('errors!');
    return;
end

press = 0;

for i = 1:n
    E = [ E0(1:i-1,:);  E0(i+1:n,:) ];
    F = [ F0(1:i-1,:);  F0(i+1:n,:) ];
  
    for j = 1:h         % h步
        w = ( E'* F )/ norm( E'*F );
        W(:, j ) = w;
        
        t = E * w;
        T(:, j ) = t;     % 主成分
        
        p = E'* t / ( norm( t)^2 );
        P(:, j ) = p;
        
        r = F'* t / ( norm( t)^2 );
        R( j ) = r;       % 系数
        
        E = E - t * p';
        F = F - t * r;
    end
    
    % 求W* 的值
    for k = 1: h
        W_h_s = eye(m);
        for j = 1:(k-1)
            W_h_s = W_h_s * ( eye(m) - W(:,j) *( P(:,j) )');
        end
        
        W_h_s = W_h_s * W(:,k);         % Wh* 的值
        
        W_s(:,k) = W_h_s;
    end
    
    F0i_s = E0(i,:) * (W_s * R');      % 求F0i的拟合值
    F0i = F0( i);
    
    % Yi_Yh_i = ( F0i - F0i_s )^2
    
    press = press + ( F0i - F0i_s )^2; % 两值相减再平方

end