%Simulate Data
%y=a + bx + cx^2
n=50;
a=25.0;
b=5.6;
c=-0.25;


x=[1:n]';
yactual=[a+b*x+c*x.^2];
subplot(2,2,1)
hold on;
plot(x,yactual,'.');

%Add random normal "noise" 
sigma=5;
e=[normrnd(0,sigma,n,1)];
y=yactual+e;
plot(x,y,'Marker','*','Color','red','LineStyle','none');
axis([0 50 0 100]);
legend('Actual','Observed',2)
hold off

%Fit the Regression Model
xmat=[[ones(1,n)]',x,x.^2];
xdim=size(xmat);
p=xdim(2);

est=inv(xmat'*xmat)*(xmat'*y);
yfit=xmat*est;
sigma_est=(y-yfit)'*(y-yfit)/(n-p);

var_est=sigma_est*inv(xmat'*xmat);

se_est=sqrt(diag(var_est));

tstat=est./se_est;

subplot(2,2,2)
hold on
plot(x,y,'Marker','*','Color','red','LineStyle','none');
plot(x,yfit,'Color','green');
axis([0 50 0 100]);
legend('Observed','Fitted',2);
hold off;

subplot(2,2,3)
hold on
plot(x,yactual,'Color','blue');
plot(x,yfit,'Color','green');
axis([0 50 0 100]);
legend('Actual','Fitted',2);
hold off;


