| 9 | from scipy.special import digamma, gamma |
| 10 | |
| 11 | def get_cost(X, K, cluster_assignments, phi, alphas, mu_means, mu_covs, a, B, orig_alphas, orig_c, orig_a, orig_B): |
| 12 | N, D = X.shape |
| 13 | total = 0 |
| 14 | ln2pi = np.log(2*np.pi) |
| 15 | |
| 16 | # calculate B inverse since we will need it |
| 17 | Binv = np.empty((K, D, D)) |
| 18 | for j in xrange(K): |
| 19 | Binv[j] = np.linalg.inv(B[j]) |
| 20 | |
| 21 | # calculate expectations first |
| 22 | Elnpi = digamma(alphas) - digamma(alphas.sum()) # E[ln(pi)] |
| 23 | Elambda = np.empty((K, D, D)) |
| 24 | Elnlambda = np.empty(K) |
| 25 | for j in xrange(K): |
| 26 | Elambda[j] = a[j]*Binv[j] |
| 27 | Elnlambda[j] = D*np.log(2) - np.log(np.linalg.det(B[j])) |
| 28 | for d in xrange(D): |
| 29 | Elnlambda[j] += digamma(a[j]/2.0 + (1 - d)/2.0) |
| 30 | |
| 31 | # now calculate the log joint likelihood |
| 32 | # Gaussian part |
| 33 | # total -= N*D*ln2pi |
| 34 | # total += 0.5*Elnlambda.sum() |
| 35 | # for j in xrange(K): |
| 36 | # # total += 0.5*Elnlambda[j] # vectorized |
| 37 | # for i in xrange(N): |
| 38 | # if cluster_assignments[i] == j: |
| 39 | # diff_ij = X[i] - mu_means[j] |
| 40 | # total -= 0.5*( diff_ij.dot(Elambda[j]).dot(diff_ij) + np.trace(Elambda[j].dot(mu_covs[j])) ) |
| 41 | |
| 42 | # mixture coefficient part |
| 43 | # total += Elnpi.sum() |
| 44 | |
| 45 | # use phi instead |
| 46 | for j in xrange(K): |
| 47 | for i in xrange(N): |
| 48 | diff_ij = X[i] - mu_means[j] |
| 49 | inside = Elnlambda[j] - D*ln2pi |
| 50 | inside += -diff_ij.dot(Elambda[j]).dot(diff_ij) - np.trace(Elambda[j].dot(mu_covs[j])) |
| 51 | # inside += Elnpi[j] |
| 52 | total += phi[i,j]*(0.5*inside + Elnpi[j]) |
| 53 | |
| 54 | |
| 55 | # E{lnp(mu)} - based on original prior |
| 56 | for j in xrange(K): |
| 57 | E_mu_dot_mu = np.trace(mu_covs[j]) + mu_means[j].dot(mu_means[j]) |
| 58 | total += -0.5*D*np.log(2*np.pi*orig_c) - 0.5*E_mu_dot_mu/orig_c |
| 59 | |
| 60 | # print "total:", total |
| 61 | |
| 62 | # E{lnp(lambda)} - based on original prior |
| 63 | for j in xrange(K): |
| 64 | total += (orig_a[j] - D - 1)/2.0*Elnlambda[j] - 0.5*np.trace(orig_B[j].dot(Elambda[j])) |
| 65 | # print "total 1:", total |
| 66 | total += -orig_a[j]*D/2.0*np.log(2) + 0.5*orig_a[j]*np.log(np.linalg.det(orig_B[j])) |
| 67 | # print "total 2:", total |
| 68 | total -= D*(D-1)/4.0*np.log(np.pi) |