update

Author: bob7783

unsupervised_class/gmm.py

@@ -25,7 +25,7 @@ def gmm(X, K, max_iter=20, smoothing=1e-2):
         M[k] = X[np.random.choice(N)]
         C[k] = np.eye(D)
 
-    costs = np.zeros(max_iter)
+    costs = []
     weighted_pdfs = np.zeros((N, K)) # we'll use these to store the PDF value of sample n and Gaussian k
     for i in range(max_iter):
         # step 1: determine assignments / resposibilities
@@ -57,7 +57,8 @@ def gmm(X, K, max_iter=20, smoothing=1e-2):
             # C[k] = np.sum(R[n,k]*np.outer(X[n] - M[k], X[n] - M[k]) for n in range(N)) / Nk + np.eye(D)*smoothing
 
 
-        costs[i] = np.log(weighted_pdfs.sum(axis=1)).sum()
+        c = np.log(weighted_pdfs.sum(axis=1)).sum()
+        costs.append(c)
         if i > 0:
             if np.abs(costs[i] - costs[i-1]) < 0.1:
                 break

unsupervised_class/kmeans.py

@@ -9,6 +9,7 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
+from sklearn.metrics.pairwise import pairwise_distances
 
 
 def d(u, v):
@@ -30,7 +31,7 @@ def cost(X, R, M):
     return cost
 
 
-def plot_k_means(X, K, max_iter=20, beta=1.0, show_plots=True):
+def plot_k_means(X, K, max_iter=20, beta=3.0, show_plots=False):
     N, D = X.shape
     M = np.zeros((K, D))
     # R = np.zeros((N, K))
@@ -40,27 +41,41 @@ def plot_k_means(X, K, max_iter=20, beta=1.0, show_plots=True):
     for k in range(K):
         M[k] = X[np.random.choice(N)]
 
-    costs = np.zeros(max_iter)
+    costs = []
+    k = 0
     for i in range(max_iter):
+        k += 1
         # step 1: determine assignments / resposibilities
         # is this inefficient?
         for k in range(K):
             for n in range(N):
-                # R[n,k] = np.exp(-beta*d(M[k], X[n])) / np.sum( np.exp(-beta*d(M[j], X[n])) for j in range(K) )
                 exponents[n,k] = np.exp(-beta*d(M[k], X[n]))
-
         R = exponents / exponents.sum(axis=1, keepdims=True)
-        # assert(np.abs(R - R2).sum() < 1e-10)
+
 
         # step 2: recalculate means
-        for k in range(K):
-            M[k] = R[:,k].dot(X) / R[:,k].sum()
+        # decent vectorization
+        # for k in range(K):
+        #     M[k] = R[:,k].dot(X) / R[:,k].sum()
+        # oldM = M
 
-        costs[i] = cost(X, R, M)
+        # full vectorization
+        M = R.T.dot(X) / R.sum(axis=0, keepdims=True).T
+        # print("diff M:", np.abs(M - oldM).sum())
+
+        c = cost(X, R, M)
+        costs.append(c)
         if i > 0:
-            if np.abs(costs[i] - costs[i-1]) < 1e-5:
+            if np.abs(costs[-1] - costs[-2]) < 1e-5:
                 break
 
+        if len(costs) > 1:
+            if costs[-1] > costs[-2]:
+                pass
+                # print("cost increased!")
+                # print("M:", M)
+                # print("R.min:", R.min(), "R.max:", R.max())
+
     if show_plots:
         plt.plot(costs)
         plt.title("Costs")
@@ -71,6 +86,7 @@ def plot_k_means(X, K, max_iter=20, beta=1.0, show_plots=True):
         plt.scatter(X[:,0], X[:,1], c=colors)
         plt.show()
 
+    print("Final cost", costs[-1])
     return M, R
 
 
@@ -98,13 +114,19 @@ def main():
     plt.show()
 
     K = 3 # luckily, we already know this
-    plot_k_means(X, K)
+    plot_k_means(X, K, beta=1.0, show_plots=True)
+
+    K = 3 # luckily, we already know this
+    plot_k_means(X, K, beta=3.0, show_plots=True)
+
+    K = 3 # luckily, we already know this
+    plot_k_means(X, K, beta=10.0, show_plots=True)
 
     K = 5 # what happens if we choose a "bad" K?
-    plot_k_means(X, K, max_iter=30)
+    plot_k_means(X, K, max_iter=30, show_plots=True)
 
     K = 5 # what happens if we change beta?
-    plot_k_means(X, K, max_iter=30, beta=0.3)
+    plot_k_means(X, K, max_iter=30, beta=0.3, show_plots=True)
 
 
 if __name__ == '__main__':

unsupervised_class/kmeans_fail.py

@@ -35,19 +35,19 @@ def donut():
 def main():
     # donut
     X = donut()
-    plot_k_means(X, 2)
+    plot_k_means(X, 2, beta=0.1, show_plots=True)
 
     # elongated clusters
     X = np.zeros((1000, 2))
     X[:500,:] = np.random.multivariate_normal([0, 0], [[1, 0], [0, 20]], 500)
     X[500:,:] = np.random.multivariate_normal([5, 0], [[1, 0], [0, 20]], 500)
-    plot_k_means(X, 2)
+    plot_k_means(X, 2, beta=0.1, show_plots=True)
 
     # different density
     X = np.zeros((1000, 2))
     X[:950,:] = np.array([0,0]) + np.random.randn(950, 2)
     X[950:,:] = np.array([3,0]) + np.random.randn(50, 2)
-    plot_k_means(X, 2)
+    plot_k_means(X, 2, show_plots=True)