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';