%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Q1.1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b = 100;
n=20;
niter = 30;
MLE = zeros(niter,1);
MAP = zeros(niter,1);
l = 0.2
for a=1:niter
    X=exprnd(1/l,n,50);
    MLE(a) = mean((1./mean(X)-l).^2);
    MAP(a) = mean(((n+a-1)./(sum(X)+b)-l).^2);
end
plot(1:niter,MLE,1:niter,MAP)
legend(['MLE', 'MAP'])

b = 100;
a=30
niter = 1000;
MLE = zeros(niter,1);
MAP = zeros(niter,1);
l = 0.2
for n=1:niter
    X=exprnd(1/l,n,50);
    MLE(n) = mean((1./mean(X)'-l).^2);
    MAP(n) = mean(((n+a-1)./(sum(X)'+b)-l).^2);
end
semilogy(1:niter,MLE,1:niter,MAP)
legend(['MLE', 'MAP'])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Q1.2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('estimate.mat')
niters = 100;
alpha = 1;
alphas = zeros(niters,1);
alpha0 = 1;
alphas0 = zeros(niters,1);
obj = zeros(niters,1);
n = size(xs,1);
for i=1:niters
    g0 = n*log(alpha0/mean(xs)) + n- n*digamma(alpha0) + sum(log(xs)) - sum(xs)/mean(xs);
    gp = n/alpha0 - n*trigamma(alpha0);
    alpha0 -= g0/gp; %Newton's method
    alphas0(i) = alpha0;

    g = n*log(alpha/mean(xs)) + n- n*digamma(alpha) + sum(log(xs)) - sum(xs)/mean(xs);
    alpha += 0.01*g; %GD
    alphas(i) = alpha;
end
plot(1:niters,alphas,1:niters,alphas0)
legend("GD","Newton")

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Q2.2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ model ] = NaiveBayes(XTrain, yTrain)
    yTrain = indicator(yTrain);
    yCounts = sum(yTrain);
    model.Py = yCounts/sum(yCounts);

    % For each p(x|y) compute mu, sigma^2 estimates, indexed by [y][x]
    model.mean = zeros(size(yTrain,2),size(XTrain,2));
    model.sigma = zeros(size(yTrain,2),size(XTrain,2));
    for i=1:size(yTrain,2)
        Xs = XTrain(logical(yTrain(:,i)),:);
        model.mean(i,:) = mean(Xs);
        model.sigma(i,:) = std(Xs);
    end
end

function [predictions] = NaiveBayesClassify(model, XTest)
log_p = zeros(size(XTest,1),size(model.mean,1));
for i=1:size(XTest,1)
  for k=1:size(model.mean,1)
    log_p(i,k) = log(model.Py(k)) + sum(log(normpdf(XTest(i,:),model.mean(k,:),model.sigma(k,:))));
  end
end
[max_values predictions] = max(log_p,[],2);
end

function [ model ] = AvgNaiveBayes(XTrain, yTrain)
    yTrain = indicator(yTrain);
    yCounts = sum(yTrain);
    model.Py = yCounts/sum(yCounts);

    model.mean = zeros(size(yTrain,2),size(XTrain,2));
    model.sigma = zeros(size(yTrain,2),size(XTrain,2));
    for i=1:size(yTrain,2)
        Xs = XTrain(logical(yTrain(:,i)),:);
        model.mean(i,:) = mean(Xs);
        model.sigma(i,:) = sqrt(var(Xs));
    end
    model.x_mean = mean(XTrain);
    model.x_sigma = std(XTrain);

    log_px = zeros(size(XTrain,2),1);
    [yMax yIndices] = max(yTrain,[],2);
    log_cond_px = zeros(size(XTrain,2),1);
    for j=1:size(XTrain,2)
      log_px(j) = sum(log(normpdf(XTrain(:,j),model.x_mean(j),model.x_sigma(j))));
      log_cond_px(j) = sum(log(normpdf(XTrain(:,j),model.mean(yIndices,j),model.sigma(yIndices,j))));
    end
    model.log_px = log_px;
    model.log_cond_px = log_cond_px;
end

function [predictions] = AvgNaiveBayesClassify_T(model, XTest, Beta)
log_p = zeros(size(XTest,1),size(model.mean,1));
for i=1:size(XTest,1)
  for k=1:size(model.mean,1)
    f = model.log_px + log(normpdf(XTest(i,:),model.x_mean,model.x_sigma))';
    g = model.log_cond_px + log(normpdf(XTest(i,:),model.mean(k,:),model.sigma(k,:)))';
    m = max(f,g);
    log_p(i,k) = log(model.Py(k)) + sum(m + log(exp(f-m) + (1/Beta)*exp(g-m)));
  end
end
[max_values predictions] = max(log_p,[],2);
end