Commented and slightly better variable naming version of posted "public" HW1, #2 solutions:

plot_lwlr.m:

function plot_lwlr(X, y, tau, res)

x = zeros(2,1);

for i=1:res,

  for j=1:res,

    x(1) = 2*(i-1)/(res-1) - 1;

    x(2) = 2*(j-1)/(res-1) - 1;

    predict(j,i) = lwlr(X, y, x, tau); %new prediction mappings

    % These loops create res-by-res "cells" between -1 to 1 (our range of x's),

    %  outputs a 0 or 1 prediction.

  end

end

figure(1);

clf;

axis off;

hold on;

colors = [-0.4 1.3];

imagesc(predict, colors);

plot((res/2)*(1+X(y==0,1))+0.5, (res/2)*(1+X(y==0,2))+0.5, 'ko');

plot((res/2)*(1+X(y==1,1))+0.5, (res/2)*(1+X(y==1,2))+0.5, 'kx');

axis equal;

axis square;

text(res/2 - res/7, res + res/20, ['tau = ' num2str(tau)], 'FontSize', 18);

---------------------------------------------------------

lwlr.m:

function y = lwlr(X_train, y_train, x, tau)

m = size(X_train,1);

n = size(X_train,2);

theta = zeros(n,1);

% compute weights, sum(_, 2) means sum the n cols of (Xi-xi)^2 into 1 col vec.

w = exp(-sum((X_train - repmat(x', m, 1)).^2, 2) / (2*tau));

% perform Newton's method

grd = ones(n,1);

lam = 1e-4;

while (norm(grd) > 1e-6)

  % hx (h(x)) equation from p.16, all others from public solution p.2 

  hx = 1 ./ (1 + exp(-X_train * theta));

  grd = X_train' * (w.*(y_train - hx)) - lam*theta;

  D = -diag(w.*hx.*(1-hx));

  H = X_train' * D * X_train -lam*eye(n);

  theta = theta - H\grd;  % same as -inv(H)*grd

end

% return predicted y

y = double(x'*theta > 0);