working kinda?

assignments/hwk03/EMG.m

@@ -22,12 +22,12 @@ function [h, m, Q] = EMG(x, k, epochs, flag)
   Q = zeros(epochs*2,1); % vector that can hold complete data log-likelihood after each E and M step
 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  % Initialise cluster means using k-means
+  % TODO: Initialise cluster means using k-means
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   [~, ~, ~, D] = kmeans(x, k);
  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  % Determine the b values for all data points
+  % TODO: Determine the b values for all data points
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   for i = 1:num_data
     row = D(i,:);
@@ -36,7 +36,7 @@ function [h, m, Q] = EMG(x, k, epochs, flag)
   end
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  % Initialize pi's (mixing coefficients)
+  % TODO: Initialize pi's (mixing coefficients)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   pi = zeros(k, 1);
   for i = 1:k
@@ -44,8 +44,8 @@ function [h, m, Q] = EMG(x, k, epochs, flag)
   end
 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  % Initialize the covariance matrix estimate
+  % TODO: Initialize the covariance matrix estimate
-  % further modifications will need to be made when doing 2(d)
+  %       further modifications will need to be made when doing 2(d)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   m = zeros(k, dim);
   for i = 1:k
@@ -63,25 +63,21 @@ function [h, m, Q] = EMG(x, k, epochs, flag)
     [h] = E_step(x, h, pi, m, S, k);
     
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    % Store the value of the complete log-likelihood function
+    % TODO: Store the value of the complete log-likelihood function
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    L = 0;
+    Q(2*n - 1) = Q_step(x, m, S, k, pi, h);
-    for i = 1:num_data
+
-      for j = 1:k
-        prior = mvnpdf(x, m(j, :), S(:, :, j));
-        L = L + h(i, j) * (log(pi(j)) + log(prior(j)));
-      end
-    end
 
     %%%%%%%%%%%%%%%%
     % M-step
     %%%%%%%%%%%%%%%%
     fprintf('M-step, epoch #%d\n', n);
-    [Q, S, m] = M_step(x, h, S, k, flag);
+    [S, m, pi] = M_step(x, h, S, k, flag);
     
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % TODO: Store the value of the complete log-likelihood function
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    Q(2*n) = Q_step(x, m, S, k, pi, h);
   end
 
 

assignments/hwk03/E_step.m

@@ -14,17 +14,17 @@ function [h] = E_step(x, h, pi, m, S, k)
   [num_data, ~] = size(x);
 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  % TODO: perform E-step of EM algorithm
+  % perform E-step of EM algorithm
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   parts = zeros(num_data, k);
 
-  for j = 1:k
+  for i = 1:k
-    parts(:, j) = pi(j) * mvnpdf(x, m(j, :), S(:, :, j));
+    parts(:, i) = pi(i) * mvnpdf(x, m(i, :), S(:, :, i));
   end
 
   s = sum(parts);
-  for i = 1:num_data
+  for j = 1:num_data
-    h(i, :) = parts(i, :) ./ s;
+    h(j, :) = parts(j, :) ./ s;
   end
 
 end

assignments/hwk03/M_step.m

@@ -20,25 +20,51 @@ function [S, m, pi] = M_step(x, h, S, k, flag)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % update mixing coefficients
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  pi = zeros(k, 1);
+  N_i = zeros(k, 1);
-  for i = 1:num_data
+  m = zeros(k, dim);
-    row = h(i, :);
+  for i = 1:k
-    maxValue = max(row);
+    N_i(i) = sum(h(:, i));
-    maxIdx = find(row == maxValue);
+    for j = 1:num_data
-    pi(maxIdx) = pi(maxIdx) + 1;
+      m(i, :) = m(i, :) + h(j, i) * x(j, :);
+    end
+  end
+  pi = N_i / num_data;
+  
+  for i = 1:k
+    m(i, :) = m(i, :) ./ N_i(i);
   end
-  pi = pi ./ num_data;
   
-
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  % TODO: update cluster means
+  % update cluster means
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
+  % m = zeros(k, dim);
+  % m = h' * x ./ N_i;
+  % for i = 1:k
+  %   m(i, :) = sum(h(:, i) .* x(i, :)) / N_i(i);
+  % end
 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  % TODO: Calculate the covariance matrix estimate
+  % Calculate the covariance matrix estimate
-  %       further modifications will need to be made when doing 2(d)
+  % further modifications will need to be made when doing 2(d)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  S = zeros(dim, dim, k) + eps;
+  for i = 1:k
+    s = zeros(dim, dim);
+    for j = 1:num_data
+      s = s + h(j, i) * (x(j, :) - m(i, :))' * (x(j, :) - m(i, :));
+    end
+    s = s / N_i(i);
+    % s = (x - m(i, :))' * ((x - m(i, :)) .* h(:, i)) / N_i(i);
+    
+    % % MAKE IT SYMMETRIC https://stackoverflow.com/a/38730499
+    % S(:, :, i) = (s + s') / 2;
+    % https://www.mathworks.com/matlabcentral/answers/366140-eig-gives-a-negative-eigenvalue-for-a-positive-semi-definite-matrix#answer_290270
+    s = (s + s') / 2;
+    % https://www.mathworks.com/matlabcentral/answers/57411-matlab-sometimes-produce-a-covariance-matrix-error-with-non-postive-semidefinite#answer_69524
+    [V, D] = eig(s);
+    s = V * max(D, eps) / V;
 
+    S(:, :, i) = s;
+  end
 
 end

assignments/hwk03/Problem2.m

@@ -26,25 +26,26 @@ function [] = Problem2()
   figure();
   for k = 4:4:12 
     fprintf("k=%d\n", k);
-
+  
     % call EM on data
     [h, m, Q] = EMG(stadium_x, k, epochs, false);
-
+  
     % get compressed version of image
     [~,class_index] = max(h,[],2);
     compress = m(class_index,:);
-
+  
     % 2(a), plot compressed image
     subplot(3,2,index)
     imagesc(permute(reshape(compress, [width, height, depth]),[2 1 3]))
     index = index + 1;
-
+  
     % 2(b), plot complete data likelihood curves
     subplot(3,2,index)
     x = 1:size(Q);
     c = repmat([1 0 0; 0 1 0], length(x)/2, 1);
     scatter(x,Q,20,c); 
     index = index + 1;
+    pause;
   end
   shg
 
@@ -69,6 +70,7 @@ function [] = Problem2()
   % TODO: plot goldy image after using clusters from k-means
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % begin code here
+  [~, ~, ~, D] = kmeans(goldy_x, k);
 
   % end code here 
   shg 

assignments/hwk03/Q_step.m

@@ -0,0 +1,11 @@
+function [LL] = Q_step(x, m, S, k, pi, h)
+  [num_data, ~] = size(x);
+  LL = 0;
+  for i = 1:k
+    N = mvnpdf(x, m(i, :), S(:, :, i));
+    for j = 1:num_data
+      LL = LL + h(j, i) * (log(pi(i) + eps) + log(N(j) + eps));
+    end
+  end
+end
+

assignments/hwk03/hw3_sol.typ

@@ -66,4 +66,16 @@ Updates:
     &= - sum_t (r^t-y^t) (v_1 frac(diff z^t_1, diff w_j) + v_2 frac(diff z^t_2, diff w_j))  \
     &= - sum_t (r^t-y^t) (x^t_j v_1 cases(0 "if" ww dot xx < 0, 1 "otherwise") + x^t_j v_2 (1 - tanh^2 (ww dot xx)))  \
     &= - sum_t (r^t-y^t) x^t_j (v_1 cases(0 "if" ww dot xx < 0, 1 "otherwise") + v_2 (1 - tanh^2 (ww dot xx)))  \
-  $
+  $
+
+#pagebreak()
+
+= Problem 2a + 2b
+
+#image("2a.png")
+
+= Problem 2c
+
+= Problem 2d
+
+MLE of $Sigma_i$