| 32 | |
| 33 | |
| 34 | def plot_k_means(X, K, max_iter=20, beta=3.0, show_plots=False): |
| 35 | N, D = X.shape |
| 36 | # R = np.zeros((N, K)) |
| 37 | exponents = np.empty((N, K)) |
| 38 | |
| 39 | # initialize M to random |
| 40 | initial_centers = np.random.choice(N, K, replace=False) |
| 41 | M = X[initial_centers] |
| 42 | |
| 43 | costs = [] |
| 44 | k = 0 |
| 45 | for i in range(max_iter): |
| 46 | k += 1 |
| 47 | # step 1: determine assignments / resposibilities |
| 48 | # is this inefficient? |
| 49 | for k in range(K): |
| 50 | for n in range(N): |
| 51 | exponents[n,k] = np.exp(-beta*d(M[k], X[n])) |
| 52 | R = exponents / exponents.sum(axis=1, keepdims=True) |
| 53 | |
| 54 | |
| 55 | # step 2: recalculate means |
| 56 | # decent vectorization |
| 57 | # for k in range(K): |
| 58 | # M[k] = R[:,k].dot(X) / R[:,k].sum() |
| 59 | # oldM = M |
| 60 | |
| 61 | # full vectorization |
| 62 | M = R.T.dot(X) / R.sum(axis=0, keepdims=True).T |
| 63 | # print("diff M:", np.abs(M - oldM).sum()) |
| 64 | |
| 65 | c = cost(X, R, M) |
| 66 | costs.append(c) |
| 67 | if i > 0: |
| 68 | if np.abs(costs[-1] - costs[-2]) < 1e-5: |
| 69 | break |
| 70 | |
| 71 | if len(costs) > 1: |
| 72 | if costs[-1] > costs[-2]: |
| 73 | pass |
| 74 | # print("cost increased!") |
| 75 | # print("M:", M) |
| 76 | # print("R.min:", R.min(), "R.max:", R.max()) |
| 77 | |
| 78 | if show_plots: |
| 79 | plt.plot(costs) |
| 80 | plt.title("Costs") |
| 81 | plt.show() |
| 82 | |
| 83 | random_colors = np.random.random((K, 3)) |
| 84 | colors = R.dot(random_colors) |
| 85 | plt.scatter(X[:,0], X[:,1], c=colors) |
| 86 | plt.show() |
| 87 | |
| 88 | print("Final cost", costs[-1]) |
| 89 | return M, R |
| 90 | |
| 91 | |