% 'hw2_prob2.m'
% 6.867 Machine Learning 
% Fall 2003

% initialization data
x = [0 0 0 0 0 1 1 1 1 1]';
y = [0 0 1 0 1 1 1 1 1 0]';

% initialize counters.
n0 = sum(x==0);
n1 = sum(x==1);
n01 = sum(x==0 & y==1);
n11 = sum(x==1 & y==1);

data = zeros(500,6);

active = 0;  % set to 0 for part (2-5)

% active learning loop...
for k = 1:500

  % calculate ML theta params
  theta0 = n01/n0;
  theta1 = n11/n1;

  % select x to sample next
  sigma0 = theta0 * (1-theta0);
  sigma1 = theta1 * (1-theta1);
  V0 = sigma0 / (n0*(n0+1));
  V1 = sigma1 / (n1*(n1+1));

  if (active)
    x = (V1 > V0);
  else
    x = (rand(1,1)>0.5); % stochastic sampling
  end

  % compute J's (for plots)
  J_est = sigma0/n0 + sigma1/n1;
  J_true = (.4*.6)/n0 + (.9*.1)/n1; % for comparison only

  % sample output for selected x
  y = query(x);    

  % update counters
  if (x) 
    n1=n1+1;
    if (y)
      n11=n11+1;
    end
  else % x==0
    n0=n0+1;
    if (y)
      n01=n01+1;
    end
  end

  % store data for generating plots
  data(k,1) = n0;
  data(k,2) = n1;
  data(k,3) = theta0;
  data(k,4) = theta1;
  data(k,5) = J_est;
  data(k,6) = J_true;

end

k = [1:500];

figure(1);
plot(k,data(k,3),'-',k,data(k,4),'--');
legend('theta0','theta1');
xlabel('iteration');
ylabel('estimate');
title('theta estimates');

figure(2);
plot(k,data(k,5),'--',k,data(k,6),'-');
legend('estimated','actual');
xlabel('iterations');
ylabel('J value');
title('overall variance');

figure(3);
n = k + 10;
plot(k,n.*data(k,5)','--',k,n.*data(k,6)','-');
legend('estimated','actual');
xlabel('iterations');
ylabel('n times J');
title('rescaled overall variance');