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)