| 96 | |
| 97 | |
| 98 | def gmm(X, K, max_iter=100): |
| 99 | N, D = X.shape |
| 100 | |
| 101 | # parameters for pi, mu, and precision |
| 102 | alphas = np.ones(K, dtype=np.float32) # prior parameter for pi (dirichlet) |
| 103 | orig_alphas = np.ones(K, dtype=np.float32) # prior parameter for pi (dirichlet) |
| 104 | # mu_means = np.zeros((K, D), dtype=np.float32) # prior mean for mu (normal) ### No! |
| 105 | # mu_covs = np.empty((K, D, D), dtype=np.float32) # prior covariance for mu (normal) |
| 106 | |
| 107 | orig_c = 10.0 |
| 108 | # for k in xrange(K): |
| 109 | # mu_covs[k] = np.eye(D)*orig_c |
| 110 | |
| 111 | |
| 112 | orig_a = np.ones(K, dtype=np.float32)*D |
| 113 | a = np.ones(K, dtype=np.float32)*D # prior for precision (wishart) |
| 114 | orig_B = np.empty((K, D, D)) |
| 115 | B = np.empty((K, D, D)) # precision (wishart) |
| 116 | empirical_cov = np.cov(X.T) |
| 117 | for k in xrange(K): |
| 118 | B[k] = (D/10.0)*empirical_cov |
| 119 | orig_B[k] = (D/10.0)*empirical_cov |
| 120 | |
| 121 | |
| 122 | # try random init instead |
| 123 | # mu_means = np.random.randn(K, D)*orig_c |
| 124 | mu_means = np.empty((K, D)) |
| 125 | for j in xrange(K): |
| 126 | mu_means[j] = X[np.random.choice(N)] |
| 127 | mu_covs = wishart.rvs(df=orig_a[0], scale=np.linalg.inv(B[0]), size=K) |
| 128 | |
| 129 | costs = np.zeros(max_iter) |
| 130 | for iter_idx in xrange(max_iter): |
| 131 | # calculate q(c[i]) |
| 132 | # phi = np.empty((N,K)) # index i = sample, index j = cluster |
| 133 | t1 = np.empty(K) |
| 134 | t2 = np.empty((N,K)) |
| 135 | t3 = np.empty(K) |
| 136 | t4 = np.empty(K) |
| 137 | |
| 138 | # calculate this first because we will use it multiple times |
| 139 | Binv = np.empty((K, D, D)) |
| 140 | for j in range(K): |
| 141 | Binv[j] = np.linalg.inv(B[j]) |
| 142 | |
| 143 | for j in xrange(K): |
| 144 | # calculate t1 |
| 145 | t1[j] = -np.log(np.linalg.det(B[j])) |
| 146 | for d in xrange(D): |
| 147 | t1[j] += digamma( (1 - d + a[j])/2.0 ) |
| 148 | |
| 149 | # calculate t2 |
| 150 | for i in xrange(N): |
| 151 | diff_ij = X[i] - mu_means[j] |
| 152 | t2[i,j] = diff_ij.dot( (a[j]*Binv[j] ).dot(diff_ij) ) |
| 153 | |
| 154 | # calculate t3 |
| 155 | t3[j] = np.trace( a[j]*Binv[j].dot(mu_covs[j]) ) |