%%%%%%%   hw1q1.m  %%%%%%%
load('q1x.dat');
load('q1y.dat');

q1x = [ones(size(q1x,1),1) q1x];
[theta, ll] = log_regression(q1x,q1y);

m=size(q1x,1);
figure; hold on;

for i=1:m
   if(q1y(i)==0)
      plot(q1x(i,2),q1x(i,3),'rx');
   else
      plot(q1x(i,2),q1x(i,3),'go');
   end
end

x = min(q1x(:,2)):.01:max(q1x(:,2));
y = -theta(1)/theta(3)-theta(2)/theta(3)*x;

plot(x,y);
xlabel('x1');
ylabel('x2');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%% log_regression.m %%%%%%%

function [theta,ll] = log_regression(X,Y)

% rows of X are training samples
% rows of Y are corresponding 0/1 values

% newton raphson: theta = theta - inv(H)* grad;
% with H = hessian, grad = gradient

m = size(X,1);
n = size(X,2);

theta = zeros(n,1);

max_iters = 50;
for i=1:max_iters
   grad = zeros(n,1);
   ll(i)=0;
   H = zeros(n,n);
   for j=1:m
      hxj = sigmoid(X(j,:)*theta);
      grad = grad + X(j,:)'*(Y(j) - hxj);
      H = H - hxj*(1-hxj)*X(j,:)'*X(j,:);
      ll(i) = ll(i) + Y(j)*log(hxj) + (1-Y(j))*log(1-hxj);
   end
   theta = theta - inv(H)*grad;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%% sigmoid.m %%%%%%%%%%%%%%
function a = sigmoid(x)

a = 1./(1+exp(-x));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
