%------------------------------------------------------------
% Tony Hyun Kim
% CS229: PS#1, Question 1b
%------------------------------------------------------------
function [theta X Y] = p1b

clear all; close all;
fSize = 16;

% Load and display data
global X Y;
X = load('q1x.dat');
Y = load('q1y.dat');
N = length(X(:,1));
X = [ones(N,1) X]; % Add a column of x0=1

sieve = (Y==1);
X0 = X(~sieve,:);
plot(X0(:,2),X0(:,3),'.'); hold on;
set(gca,'FontSize',fSize);

X1 = X(sieve,:);
plot(X1(:,2),X1(:,3),'rx');
xlabel('x_1');
ylabel('x_2');
title('PS#1, Q1, Parts b,c');
legend('y^{(i)}=0','y^{(i)}=1','Location','Best');

% Newton's method to iteratively find MLE
%   The update equation comes from pg. 21 of the course notes
Niter = 10;
theta = [0 0 0]';
for i = 1:Niter
    g = gradient(theta);
    H = hessian(theta);
    theta = theta - H\g;
end

% Draw the classification boundary where h(x)=0.5
%   Based on the defn of sigmoid(z), it is clear that the 
%   desired boundary lies where z=0.
x1 = linspace(0,8);
x2 = -1/theta(3)*(theta(1)+theta(2)*x1);
plot(x1,x2,'k');

%------------------------------------------------------------
% Helper functions
%   Matrix representations for the gradient and hessian of
%   the log likelihood follows my notation from part (a)
%------------------------------------------------------------
function g = sigmoid(z)
g = 1./(1+exp(-z));

% The following function returns the gradient of
%   the log likelihood l(theta), where 
%   theta is a (n+1)-dim column vector
function g = gradient(theta)
global X Y;
e = Y - sigmoid(X*theta);
g = X'*e;

% The following function returns the Hessian of
%    the log likelihood l(theta) as a 
%    function of theta
function H = hessian(theta)
global X;
h = sigmoid(X*theta);
F = diag(h.*(1-h));
H = -X'*F*X;