| 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): |