| 77 | |
| 78 | |
| 79 | def run(num=1, T=500): |
| 80 | X = pd.read_csv('X_set%s.csv' % num, header=None).as_matrix() |
| 81 | Y = pd.read_csv('y_set%s.csv' % num, header=None).as_matrix().flatten() |
| 82 | Z = pd.read_csv('z_set%s.csv' % num, header=None).as_matrix().flatten() |
| 83 | N, D = X.shape |
| 84 | print X.shape, Y.shape, Z.shape |
| 85 | |
| 86 | a0 = 1e-16 |
| 87 | b0 = 1e-16 |
| 88 | e0 = 1 |
| 89 | f0 = 1 |
| 90 | |
| 91 | # params for q(w) - doesn't matter what we set it to, we'll update this first |
| 92 | C = np.eye(D) |
| 93 | mu = np.zeros(D) |
| 94 | |
| 95 | # params for q(lambda) |
| 96 | e = e0 |
| 97 | f = f0 |
| 98 | |
| 99 | # params for q(alpha) |
| 100 | a = np.ones(D)*a0 |
| 101 | b = np.ones(D)*b0 |
| 102 | a0ones = np.ones(D)*a0 |
| 103 | |
| 104 | # objective |
| 105 | L = np.empty(T) |
| 106 | |
| 107 | for t in xrange(T): |
| 108 | # update q(w) |
| 109 | C = np.linalg.inv(np.diag(1.0*a/b) + (1.0*e/f)*X.T.dot(X)) |
| 110 | mu = C.dot((1.0*e/f)*X.T.dot(Y)) |
| 111 | |
| 112 | # update q(alpha) |
| 113 | a = a0ones + 0.5 |
| 114 | b = b0 + 0.5*(np.diag(C) + mu*mu) |
| 115 | # for k in xrange(D): |
| 116 | # a[k] = a0 + 0.5 |
| 117 | # b[k] = b0 + 0.5*(C[k,k] + mu[k]*mu[k]) |
| 118 | |
| 119 | # update q(lambda) |
| 120 | e = e0 + N/2.0 |
| 121 | sum_for_f = 0 |
| 122 | # for i in xrange(N): |
| 123 | # delta = Y[i] - X[i].dot(mu) |
| 124 | # sum_for_f += delta*delta + X[i].dot(C).dot(X[i]) |
| 125 | delta = Y - X.dot(mu) |
| 126 | sum_for_f = delta.dot(delta) + np.trace(X.dot(C).dot(X.T)) |
| 127 | f = f0 + 0.5*sum_for_f |
| 128 | |
| 129 | # update L |
| 130 | L[t] = objective(X, Y, C, mu, a, b, e, f, a0, b0, e0, f0) |
| 131 | if t % 20 == 0: |
| 132 | print "t:", t |
| 133 | if num == 3: |
| 134 | print "L:", L[t] |
| 135 | |
| 136 | # plot 1/E[alpha] |