| 15 | return np.log(b) - a - np.log(np.abs(gamma(a))) + (a - 1)*digamma(a) |
| 16 | |
| 17 | def objective(X, Y, C, mu, a, b, e, f, a0, b0, e0, f0): |
| 18 | log2pi = np.log(2*np.pi) |
| 19 | N, D = X.shape |
| 20 | |
| 21 | # E(lnX) = digamma(a) - ln(b) for X ~ Gamma(a,b) |
| 22 | E_ln_lambda = digamma(e) - np.log(f) |
| 23 | E_ln_alpha = digamma(a) - np.log(b) |
| 24 | |
| 25 | # model likelihood |
| 26 | total = (N/2.0)*(E_ln_lambda - log2pi) |
| 27 | data_total = 0 |
| 28 | for i in xrange(N): |
| 29 | delta = Y[i] - X[i].dot(mu) |
| 30 | data_total += delta*delta + X[i].dot(C).dot(X[i]) |
| 31 | total -= (float(e)/f)/2.0 * data_total |
| 32 | |
| 33 | # print "total after model likelihood:", total |
| 34 | |
| 35 | # w likelihood |
| 36 | total -= (D/2.0)*log2pi |
| 37 | for k in xrange(D): |
| 38 | total += 0.5*(E_ln_alpha[k] - (float(a[k])/b[k])*(C[k,k] + mu[k]*mu[k])) |
| 39 | |
| 40 | # print "total after w likelihood:", total |
| 41 | |
| 42 | # lambda likelihood |
| 43 | total += e0*np.log(f0) - np.log(gamma(e0)) + (e0 - 1)*E_ln_lambda - f0*(float(e)/f) |
| 44 | |
| 45 | # print "total after lambda likelihood:", total |
| 46 | |
| 47 | # alpha likelihood |
| 48 | for k in xrange(D): |
| 49 | total += a0*np.log(b0) - np.log(gamma(a0)) + (a0 - 1)*E_ln_alpha[k] - b0*(float(a[k])/b[k]) |
| 50 | |
| 51 | # print "total after alpha likelihood:", total |
| 52 | |
| 53 | # entropy |
| 54 | # TODO: calculate this manually |
| 55 | # total -= mvn.entropy(mean=mu, cov=C) |
| 56 | # e1 = mvn.entropy(cov=C) |
| 57 | # e2 = 0.5*np.log( np.linalg.det(2*np.pi*np.e*C) ) |
| 58 | # print "e1:", e1, "e2:", e2 |
| 59 | # total += 0.5*np.log( np.linalg.det(2*np.pi*np.e*C) ) |
| 60 | |
| 61 | total += mvn.entropy(cov=C) |
| 62 | # print "det(C):", np.linalg.det(C) |
| 63 | # print "total after lnq(w):", total |
| 64 | |
| 65 | # total -= gamma_dist.entropy(e, scale=1.0/f) |
| 66 | # e3 = gamma_dist.entropy(e, scale=1.0/f) |
| 67 | # e4 = -e_ln_q_gamma(e, f) |
| 68 | # print "e3:", e3, "e4:", e4 |
| 69 | # assert(np.abs(e3 - e4) < 1e-8) |
| 70 | total += gamma_dist.entropy(e, scale=1.0/f) |
| 71 | # total -= e_ln_q_gamma(e, f) |
| 72 | # print "total after lnq(lambda):", total |
| 73 | for k in xrange(D): |
| 74 | # total -= e_ln_q_gamma(a[k], b[k]) |