| 9 | |
| 10 | |
| 11 | class ProbitRegression: |
| 12 | def fit(self, X, Y, sigma=1.5, lam=1, show_w=set(), Q=None): |
| 13 | |
| 14 | # setup |
| 15 | N, D = X.shape |
| 16 | self.w = np.random.randn(D) / np.sqrt(D) # does not work if you don't scale first! |
| 17 | Eq = np.zeros(N) |
| 18 | idx1 = (Y == 1) |
| 19 | idx0 = (Y == 0) |
| 20 | A = lam*np.eye(D) + X.T.dot(X) / sigma**2 |
| 21 | costs = [] |
| 22 | |
| 23 | for t in xrange(100): |
| 24 | # calculate ln(p(Y, w | X)) |
| 25 | cdf = norm.cdf(X.dot(self.w) / sigma) |
| 26 | cost = -lam/2 * self.w.dot(self.w) + Y.dot(np.log(cdf)) + (1 - Y).dot(np.log(1 - cdf)) |
| 27 | costs.append(cost) |
| 28 | |
| 29 | # E step |
| 30 | Xw = X.dot(self.w) |
| 31 | pdf = norm.pdf(-Xw / sigma) |
| 32 | cdf = norm.cdf(-Xw / sigma) |
| 33 | Eq[idx1] = Xw[idx1] + sigma*(pdf[idx1] / (1 - cdf[idx1])) |
| 34 | Eq[idx0] = Xw[idx0] + sigma*(-pdf[idx0] / cdf[idx0]) |
| 35 | |
| 36 | # M step |
| 37 | b = X.T.dot(Eq) / sigma**2 |
| 38 | self.w = np.linalg.solve(A, b) |
| 39 | |
| 40 | if show_w and t in show_w: |
| 41 | plot_image(self.w, Q, "iteration: %s" % (t+1)) |
| 42 | |
| 43 | |
| 44 | plt.plot(costs) |
| 45 | plt.show() |
| 46 | |
| 47 | self.sigma = sigma |
| 48 | self.lam = lam |
| 49 | |
| 50 | |
| 51 | def predict_proba(self, X): |
| 52 | N, D = X.shape |
| 53 | return norm.cdf(X.dot(self.w) / self.sigma) |
| 54 | |
| 55 | def predict(self, X): |
| 56 | return np.round(self.predict_proba(X)) |
| 57 | |
| 58 | def score(self, X, Y): |
| 59 | return np.mean(self.predict(X) == Y) |
| 60 | |
| 61 | def confusion_matrix(self, X, Y): |
| 62 | P = self.predict(X) |
| 63 | M = np.zeros((2, 2)) |
| 64 | M[0,0] = np.sum(P[Y == 0] == Y[Y == 0]) |
| 65 | M[0,1] = np.sum(P[Y == 0] != Y[Y == 0]) |
| 66 | M[1,0] = np.sum(P[Y == 1] != Y[Y == 1]) |
| 67 | M[1,1] = np.sum(P[Y == 1] == Y[Y == 1]) |
| 68 | return M |