function [Bhat VEM] = OLSregress(X,y)
% 在观测矩阵X中不必输入常数项1
% 程序会自动给出带有常数项的估计系数Bhat
[n p] = size(X);
dat = [X,y];
mdat = mean([X,y]); % 均值,用于归零和确定常数项系数
dat = dat - ones(n,1)*mdat;%归零
[U S V] = svd(dat);
A=dat(:,1:p);b=dat(:,p+1);d=diag(S); B = pinv(A'*A-d(p+1)^2*eye(p))*A'*b;Bhat = [mdat*[-B;1];B];%正交最小二乘估计系数
VEM = sum(d(1:p))/sum(d);% 线性模型能解释的方差(VEM)