%PROGRAM NAME: learning.m
% X1= Number of problems: Sx1=1:10
% X2 = Hrs of coordinated physical activity/week: Sx2=0:.5:10
% Y1 = Exam 1 score: Sy1=0:100
% Y2 = Exam 2 score: Sy2=0:100 one month later
%-------------------------------------------------------------
% X1: Uniform(a1,b1) where a1=1 & a2=10
a1=1; b1=10; mX1=(a1+b1)/2; sX1=sqrt(((b1-a1+1)^2-1)/12);
SX1=1:10; SX1=SX1'; %X1 Sample Space
% X2: Uniform(a2,b2) where a2=0 & 22=10 in 1/2 hr increments
a2=0.5; b2=20; mX2=(a2+b2)/2; sX2=sqrt(((b2-a2+1)^2-1)/12);
mX2=mX2/2; sX2=sX2/2; %Scale increments to hours
SX2=0.5:.5:10; SX2=SX2';  %X2 Sample Space
mX = [mX1 ; mX2];
cXX = [sX1^2 , 0 ; 0 sX2^2]; %X covariance matrix
%--------------------------------------------------------------
%The model for [Y1;Y2];
A = [3.5 , 1.8 ; -1.5 , 1]; b = [17.5 ; 56.5];
mY = A*mX + b;
cXY=cXX*A; %XY covariance matrix
sE1=3; sE2=6; %Model error std. deviations
cEE = [sE1^2 , 0 ; 0 , sE2^2]; % E covariance matrix
cYY=A*cXX*A' + cEE; %YY covariance matrix
% VALIDATE:
%Ahat=cXX^-1 * cXY
%bhat = mY - Ahat*mX
%THE ABOVE IS ACCURATE :)
%============================================================
%Simulate Data:
m = 10000;   % Sample Size
%Error Data:
e = mvnrnd([0,0],cEE,m);
%X Data:
nSX1=length(SX1);
rx1=fix(m/nSX1);
x1=[];
for j=1:rx1
x1j=randperm(nSX1)';
x1=[x1;x1j];
end
%----------------------
%Enable the following to see the influence of not randomizing x1:
%x1=kron(ones(rx1,1),SX1);
%----------------------
nSX2=length(SX2);
rx2=fix(m/nSX2);
x2=[];
for j=1:rx2
x2j =0.5*randperm(nSX2)';
x2=[x2;x2j];
end
%----------------------
%Enable the following to see the influence of not randomizing x2:
%x2=kron(ones(rx2,1),SX2);
%----------------------
x=[x1 , x2];
%Y Data:
y=zeros(m,2);
for k=1:m
    y(k,:)=  x(k,:)*A' + b' + e(k,:);
end
xy = [x , y];
%==========================================
mhat=mean(xy);
chat=cov(xy);
cXXhat=chat(1:2,1:2);
cXYhat=chat(3:4,1:2);
Ahat=cXYhat*cXXhat^-1;
bhat=mhat(3:4)' - Ahat*mhat(1:2)';
%----------------------------------------
%Compare to the TRUTH:
AAhat=[A , Ahat]
bbhat=[b' , bhat']
pause
cXXcXXhat=[cXX , cXXhat]
cXYcXYhat=[cXY , cXYhat]




