x=[2,2,3,4,4,5,6,7,8,9]';
y=[10,11,12,13,14,15,20,18,22,25]';
n=length(x);
m=1;
V=ones(n,2); %initalize V matrix
V(:,2)=x;
A=V'*V;
b=V'*y;
w=A\b; %solve Normal Equations
disp(w)
fit=V*w;
plot(x,fit)
hold on
plot(x,y,'*')
hold off
%standard error of estimataion and coefficient of determination
r=fit-y;
sse=r'*r;
StdErrEst=sqrt(sse/(n-m-1));
ybar=sum(y)/n;
var=y-ybar;
sst=var'*var;
ssr=sst-sse;
rsq=ssr/sst;
disp([ybar,rsq,StdErrEst])
%standard residuals
xbar=sum(x)/n;
varx=x-xbar;
sstx=varx'*varx;
for i=1:n
    h=(x(i)-xbar)^2/sstx+1/n;
    StdErrRes=StdErrEst*sqrt(1-h);
    StdRes(i)=r(i)/StdErrRes;
end
%disp(StdRes')
%inference
b1=w(2);
sumxsq=x'*x;
sb1=StdErrEst/sqrt(sumxsq-sum(x)^2/n);
tstat=b1/sb1;
pptdist=2.306; % n-2=8 95%
db1=pptdist*sb1;
b1p=b1+db1;
b1m=b1-db1;
disp([sb1,tstat,b1m,b1p])
xp=5;
yp=w(1)+xp*w(2);
hxp=1/n;
dypbar=StdErrEst*pptdist*sqrt(hxp);
dypran=StdErrEst*pptdist*sqrt(1+hxp);
ypbarp=yp+dypbar;
ypbarm=yp-dypbar;
ypranp=yp+dypran;
ypranm=yp-dypran;
disp([xp,yp,ypbarm,ypbarp,ypranm,ypranp])