% 'hw2_prob3.m'
% 6.867 Machine Learning - Fall 03
% HW2 Problem 3

% (3-1)
load data.mat;
n_train = size(trainD.y,1);
n_test = size(testD.y,1);

figure(1); plotdata(trainD); title('training data');
figure(2); plotdata(testD); title('test data');
mix = mixtrain(trainD.X,trainD.y);

% part a
% classify training data
train_loglik0 = mixloglik(trainD.X,zeros(size(trainD.y)),mix);
train_loglik1 = mixloglik(trainD.X,ones(size(trainD.y)),mix);
mix_train_err = mean(trainD.y ~= (train_loglik1 > train_loglik0))
% classify test data
test_loglik0 = mixloglik(testD.X,zeros(size(testD.y)),mix);
test_loglik1 = mixloglik(testD.X,ones(size(testD.y)),mix);
mix_test_err = mean(testD.y ~= (test_loglik1 > test_loglik0))
% part b
figure(3); mixboundary(mix,trainD); title('mix decision boundary');
% part c : compute the average conditional log-likelihoods for 
% comparison with logistic models
mix_condloglik_train = (sum(mixloglik(trainD.X,trainD.y,mix)- ...
			 log(exp(train_loglik0)+exp(train_loglik1))))/n_train
mix_condloglik_test = (sum(mixloglik(testD.X,testD.y,mix)- ...
			 log(exp(test_loglik0)+exp(test_loglik1))))/n_test

% (3-2)
% part a
w = logisticreg(trainD.X,trainD.y)
input_train = trainD.X * w(2:length(w)) + w(1); % discriminant
input_test = testD.X * w(2:length(w)) + w(1);
logistic_train_err = mean(trainD.y ~= (input_train > 0))
logistic_test_err = mean(testD.y ~= (input_test > 0))

condloglik_train = (sum(loglogistic(input_train(find(trainD.y==1)))) + ...
    sum(loglogistic(-input_train(find(trainD.y==0)))))/n_train
% where we have used the fact that log(1-g(z)) = log g(-z)

condloglik_test = (sum(loglogistic(input_test(find(testD.y==1)))) + ...
    sum(loglogistic(-input_test(find(testD.y==0)))))/n_test

% part b
figure(4); boundary(w,trainD); title('1st-order logistic decision boundary');

% (3-3)
phi = degexpand(trainD.X,2);
w = logisticreg(phi,trainD.y)
input_train = phi * w(2:length(w)) + w(1);
phi = degexpand(testD.X,2);
input_test = phi * w(2:length(w)) + w(1);
logistic_train_err2 = mean(trainD.y ~= (input_train > 0))
logistic_test_err2 = mean(testD.y ~= (input_test > 0))

condloglik_train2 = (sum(loglogistic(input_train(find(trainD.y==1)))) + ...
    sum(loglogistic(-input_train(find(trainD.y==0)))))/n_train
condloglik_test2 = (sum(loglogistic(input_test(find(testD.y==1)))) + ...
    sum(loglogistic(-input_test(find(testD.y==0)))))/n_test

figure(5); boundary(w,trainD,2); title('2nd-order logistic decision boundary');

% (3-6)
for k = 1:4
  phi = degexpand(trainD.X,k);
  w = logisticreg(phi,trainD.y);
  input_train = phi * w(2:length(w)) + w(1);

  condloglik_train = (sum(loglogistic(input_train(find(trainD.y==1)))) + ...
      sum(loglogistic(-input_train(find(trainD.y==0)))))/n_train;

  phi = degexpand(testD.X,k);
  input_test = phi * w(2:length(w)) + w(1);

  condloglik_test = (sum(loglogistic(input_test(find(testD.y==1)))) + ...
    sum(loglogistic(-input_test(find(testD.y==0)))))/n_test;

  fprintf('deg: %d, ave train log-lik: %f, ave test log-lik: %f\n',k, ...
	  condloglik_train,condloglik_test);
end