pca_regression.m

基于主元分析的多元线性回归程序,输出回归系数

function [b,pctExplained,ess]=pca_regression(x,y,k)
%基于主元分析的线性回归
%x：自变量
%y：因变量
%k：主元个数 the number of princomponent
%b：回归系数，y=b*x
%pctExplained：各主成分解释的数据变化the amount of explained variance for each component
%ess：拟合残差
%作者：秦强 2009-3-20
%参考：《偏最小二乘回归方法及其应用》，princomp命令和matlab官方网站

%例子：
%x1=3*rand(50,1);x2=3*rand(50,1);x3=2*x1-x2;x=[x1,x2,x3,ones(50,1)];y=x1+20*x2+20*x3+2;
%[b,pctExplained,ess]=pca_regression(x,y,3)
%
if nargout < 2
    k=2;
end

[u,s,v]=svd(x);

t=u(:,1:size(s,2))*diag(diag(s));%s的对角元素乘u的列
p=v;
[n,m]=size(x);

tk=t(:,1:k);%主元
pk=p(:,1:k);%负荷向量
%X=t1*p1'+t2*p2'...+tm*pm'

b=pk*(inv(tk'*tk))*(tk')*y;%回归系数

sigma = diag(s);
sigma = sigma ./ sqrt(n-1);
latent = sigma .^2;%特征根

pctExplained = latent' ./ sum(latent);%各主成分解释的数据变化

yp=x*b;
e=y-yp;
ess=norm(e);%拟合误差的平方和开根号
figure(1);
plot(e,'r-*');
b=b';