From bc30320eeff38f6178ab765d3b26eabe002e9c91 Mon Sep 17 00:00:00 2001
From: Michael Zhang <mail@mzhang.io>
Date: Thu, 16 Nov 2023 02:46:02 -0600
Subject: [PATCH] ok kinda works?

---
 assignments/hwk03/EMG.m    | 26 +++++++++++++-------------
 assignments/hwk03/E_step.m | 10 +++++-----
 assignments/hwk03/M_step.m | 38 ++++++++++++++++++++++++++------------
 3 files changed, 44 insertions(+), 30 deletions(-)

diff --git a/assignments/hwk03/EMG.m b/assignments/hwk03/EMG.m
index 93a1749..5f42882 100644
--- a/assignments/hwk03/EMG.m
+++ b/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
-  % further modifications will need to be made when doing 2(d)
+  % TODO: Initialize the covariance matrix estimate
+  %       further modifications will need to be made when doing 2(d)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   m = zeros(k, dim);
   for i = 1:k
@@ -63,21 +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;
-    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
+    % 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
 
     %%%%%%%%%%%%%%%%
     % 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
diff --git a/assignments/hwk03/E_step.m b/assignments/hwk03/E_step.m
index a7a6198..b64f19e 100644
--- a/assignments/hwk03/E_step.m
+++ b/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
-    parts(:, j) = pi(j) * mvnpdf(x, m(j, :), S(:, :, j));
+  for i = 1:k
+    parts(:, i) = pi(i) * mvnpdf(x, m(i, :), S(:, :, i));
   end
 
   s = sum(parts);
-  for i = 1:num_data
-    h(i, :) = parts(i, :) ./ s;
+  for j = 1:num_data
+    h(j, :) = parts(j, :) ./ s;
   end
 
 end
\ No newline at end of file
diff --git a/assignments/hwk03/M_step.m b/assignments/hwk03/M_step.m
index e205eec..013d037 100644
--- a/assignments/hwk03/M_step.m
+++ b/assignments/hwk03/M_step.m
@@ -21,24 +21,38 @@ function [S, m, pi] = M_step(x, h, S, k, flag)
   % update mixing coefficients
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   pi = zeros(k, 1);
-  for i = 1:num_data
-    row = h(i, :);
-    maxValue = max(row);
-    maxIdx = find(row == maxValue);
-    pi(maxIdx) = pi(maxIdx) + 1;
+  N_i = zeros(k, 1);
+  for i = 1:k
+    N_i(i) = sum(h(:, i));
   end
-  pi = pi ./ num_data;
+  pi = N_i / 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
-  %       further modifications will need to be made when doing 2(d)
+  % Calculate the covariance matrix estimate
+  % further modifications will need to be made when doing 2(d)
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
+  S = zeros(dim, dim, k);
+  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);
+    
+    % % 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;
+    [V, D] = eig(s);
+    S(:, :, i) = V * max(D,eps) / V;
+  end
 
 end
\ No newline at end of file