(X, T=500)
| 73 | |
| 74 | |
| 75 | def gmm(X, T=500): |
| 76 | N, D = X.shape |
| 77 | |
| 78 | m0 = X.mean(axis=0) |
| 79 | c0 = 0.1 |
| 80 | a0 = float(D) |
| 81 | B0 = c0*D*np.cov(X.T) |
| 82 | alpha0 = 1.0 |
| 83 | |
| 84 | # cluster assignments - originally everything is assigned to cluster 0 |
| 85 | C = np.zeros(N) |
| 86 | |
| 87 | # keep as many as we need for each gaussian |
| 88 | # originally we sample from the prior |
| 89 | # TODO: just use the function above |
| 90 | precision0 = wishart.rvs(df=a0, scale=np.linalg.inv(B0)) |
| 91 | covariances = [np.linalg.inv(precision0)] |
| 92 | means = [mvn.rvs(mean=m0, cov=covariances[0]/c0)] |
| 93 | |
| 94 | cluster_counts = [1] |
| 95 | K = 1 |
| 96 | observations_per_cluster = np.zeros((T, 6)) |
| 97 | for t in xrange(T): |
| 98 | if t % 20 == 0: |
| 99 | print t |
| 100 | # 1) calculate phi[i,j] |
| 101 | # Notes: |
| 102 | # MANY new clusters can be made each iteration |
| 103 | # A cluster can be DESTROYED if a x[i] is the only pt in cluster j and gets assigned to a new cluster |
| 104 | # phi = np.empty((N, K)) |
| 105 | list_of_cluster_indices = range(K) |
| 106 | next_cluster_index = K |
| 107 | # phi = [] # TODO: do we need this at all? |
| 108 | for i in xrange(N): |
| 109 | phi_i = {} |
| 110 | for j in list_of_cluster_indices: |
| 111 | # don't loop through xrange(K) because clusters can be created or destroyed as we loop through i |
| 112 | nj_noti = np.sum(C[:i] == j) + np.sum(C[i+1:] == j) |
| 113 | if nj_noti > 0: |
| 114 | # existing cluster |
| 115 | # phi[i,j] = N(x[i] | mu[j], cov[j]) * nj_noti / (alpha0 + N - 1) |
| 116 | # using the sampled mu / covs |
| 117 | phi_i[j] = mvn.pdf(X[i], mean=means[j], cov=covariances[j]) * nj_noti / (alpha0 + N - 1.0) |
| 118 | |
| 119 | # new cluster |
| 120 | # create a possible new cluster for every sample i |
| 121 | # but only keep it if sample i occupies this new cluster j' |
| 122 | # i.e. if C[i] = j' when we sample C[i] |
| 123 | # phi[i,j'] = alpha0 / (alpha0 + N - 1) * p(x[i]) |
| 124 | # p(x[i]) is a marginal integrated over mu and precision |
| 125 | phi_i[next_cluster_index] = alpha0 / (alpha0 + N - 1.0) * marginal(X[i], c0, m0, a0, B0) |
| 126 | |
| 127 | # normalize phi[i] and assign C[i] to its new cluster by sampling from phi[i] |
| 128 | normalize_phi_hat(phi_i) |
| 129 | |
| 130 | # if C[i] = j' (new cluster), generate mu[j'] and cov[j'] |
| 131 | C[i] = sample_cluster_identity(phi_i) |
| 132 | if C[i] == next_cluster_index: |
no test coverage detected