working kinda?

assignments/hwk03/EMG.m

@@ -65,13 +65,8 @@ function [h, m, Q] = EMG(x, k, epochs, flag)
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % TODO: Store the value of the complete log-likelihood function
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    L = 0;
-    % for i = 1:num_data
-    %   for j = 1:k
-    %     prior = mvnpdf(x, m(j, :), S(:, :, j));
-    %     L = L + h(i, j) * (log(pi(i)) + log(prior(i)));
-    %   end
-    % end
+    Q(2*n - 1) = Q_step(x, m, S, k, pi, h);
+
 
     %%%%%%%%%%%%%%%%
     % M-step
@@ -82,6 +77,7 @@ function [h, m, Q] = EMG(x, k, epochs, 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/M_step.m

@@ -20,18 +20,25 @@ function [S, m, pi] = M_step(x, h, S, k, flag)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % update mixing coefficients
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  pi = zeros(k, 1);
   N_i = zeros(k, 1);
+  m = zeros(k, dim);
   for i = 1:k
     N_i(i) = sum(h(:, i));
+    for j = 1:num_data
+      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
+  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % update cluster means
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  m = zeros(k, dim);
-  m = h' * x ./ N_i;
+  % m = zeros(k, dim);
+  % m = h' * x ./ N_i;
   % for i = 1:k
   %   m(i, :) = sum(h(:, i) .* x(i, :)) / N_i(i);
   % end
@@ -40,19 +47,24 @@ function [S, m, pi] = M_step(x, h, S, k, flag)
   % Calculate the covariance matrix estimate
   % further modifications will need to be made when doing 2(d)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-  S = zeros(dim, dim, k);
+  S = zeros(dim, dim, k) + eps;
   for i = 1:k
-    % for j = 1:num_data
-    %   S(:, :, i) = S(:, :, i) + h(j, i) * (x(j, :) - m(i, :)) * (x(j, :) - m(i, :))';
-    % end
-    s = (x - m(i, :))' * ((x - m(i, :)) .* h(:, i)) / N_i(i);
+    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(:, :, i) = V * max(D,eps) / V;
+    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$