%Program Name: humpday.m
%Truth Model:
muX=12; sigX=2;
muY=50; sigY=10;
covXY=sigX^2;
n=50;
a=covXY/sigX^2;
b=muY-a*muX;
MXY=[muX , muY];
CXY=[sigX^2 covXY ; covXY sigY^2];
%=====================================
nsim=10e2;
ahat=zeros(nsim,1);
bhat=ahat;
for k=1:nsim
    xy=mvnrnd(MXY,CXY,n);
    MXYhat=mean(xy);
    CXYhat=cov(xy);
    ahat(k)=CXYhat(1,2)/CXYhat(1,1);
    bhat(k)=MXYhat(2)-ahat(k)*MXYhat(1);
end
figure(1)
hist(ahat,20)
title('Histogram of ahat')
pause
figure(2)
hist(bhat,20)
title('Histogram of bhat')
pause
figure(3)
plot(ahat,bhat,'*')
title('Scatter Plot of ahat vs. bhat')