clear;
n=1000;
x1=randn(n,1);
x2=randn(n,1);
y1=1*(x1+2*x2)+2*randn(n,1);
y2=0.5*(x1+2*x2)+2*randn(n,1);
y3=0.05*(x1+2*x2)+2*randn(n,1);
phi_0=[0 0 1];

%Minimization
y=y1;
[phi,fval] = fminunc(@ll,phi_0,[],y,x1,x2); 
est_1=phi;
y=y2;
[phi,fval] = fminunc(@ll,phi_0,[],y,x1,x2); 
est_2=phi;
y=y3;
[phi,fval] = fminunc(@ll,phi_0,[],y,x1,x2); 
est_3=phi;


%Plots
z1=zeros(30,30);
z2=zeros(30,30);
z3=zeros(30,30);
theta=0:4/29:4;
sigmasq=1:7/29:8;
for j=1:30
    for i=1:30
z1(i,j)=-n*log(2*pi)/2 - n*log(sigmasq(j))/2 -((y1-1*(x1+theta(i)*x2))'*(y1-1*(x1+theta(i)*x2))/(2*sigmasq(j)));
    end
end
for j=1:30
    for i=1:30
z2(i,j)=-n*log(2*pi)/2 - n*log(sigmasq(j))/2 -((y2-0.5*(x1+theta(i)*x2))'*(y2-0.5*(x1+theta(i)*x2))/(2*sigmasq(j)));
    end
end
for j=1:30
    for i=1:30
z3(i,j)=-n*log(2*pi)/2 - n*log(sigmasq(j))/2 -((y3-0.05*(x1+theta(i)*x2))'*(y3-0.05*(x1+theta(i)*x2))/(2*sigmasq(j)));
    end
end

figure;
mesh(sigmasq,theta,z1);
figure;
mesh(sigmasq,theta,z2);
figure;
mesh(sigmasq,theta,z3);