ml-class / mlclass / logistic_regression.py

import warnings

from scipy.optimize import fmin_l_bfgs_b
import numpy as np
import scipy.sparse

def append_ones(X):
    m, n = X.shape
    if scipy.sparse.issparse(X):
        return scipy.sparse.hstack((np.ones((m, 1)), X)).tocsr()
    else:
        return np.column_stack((np.ones(m), X))

def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))

class LogisticRegressionGD:
    '''Gradient descent'''

    def __init__(self, lambda_=1, **fmin_args):
        self.lambda_ = lambda_
        self.fmin_args = fmin_args

    def cost_grad(self, theta, X, y):
        m = X.shape[0]
        sx = sigmoid(X.dot(theta))

        j = -np.log(np.where(y, sx, 1 - sx)).sum()
        j += self.lambda_ * theta.dot(theta) / 2.0
        j /= m

        grad = X.T.dot(sx - y) / m
        grad += self.lambda_ * theta / m

        return j, grad

    def fit(self, X, y):
        theta = np.zeros(X.shape[1])
        self.theta, j, ret = fmin_l_bfgs_b(self.cost_grad, theta, args=(X, y),
            **self.fmin_args)
        if ret['warnflag'] != 0:
            warnings.warn('L-BFGS failed to converge')

    def predict(self, X):
        return sigmoid(X.dot(self.theta))


class LogisticRegressionNR:
    '''Newton-Raphson method'''
    def __init__(self, eps=1e-7):
        self.eps = eps

    def fit(self, X, y):
        self.theta = np.zeros(X.shape[1])

        j_prev = 0
        while True:
            # compute the cost, gradient, and the Hessian
            sx = sigmoid(X.dot(self.theta))
            j = y.dot(np.log(sx)) + (1 - y).dot(np.log(1 - sx))
            grad = X.T.dot(y - sx)
            hess = -X.T.dot((sx * (1 - sx))[:,None] * X)

            if np.abs(j - j_prev) < self.eps:
                break
            self.theta -= np.linalg.inv(hess).dot(grad)
            j_prev = j

    def predict(self, X):
        return sigmoid(X.dot(self.theta))

class SoftmaxRegressionGD:
    def __init__(self, lambda_=1, **fmin_args):
        self.lambda_ = lambda_
        self.fmin_args = fmin_args

    def cost_grad(self, Theta_flat):
        m, n = self.X.shape
        k = self.Y.shape[1]

        Theta = Theta_flat.reshape((k, n))

        P = np.exp(self.X.dot(Theta.T))
        P /= np.sum(P, axis=1)[:,None]

        j = -np.sum(np.log(P) * self.Y)
        j += self.lambda_ * Theta_flat.dot(Theta_flat) / 2.0
        j /= m

        grad = self.X.T.dot(P - self.Y).T

        grad += self.lambda_ * Theta
        grad /= m

        return j, grad.ravel()

    def fit(self, X, Y):
        self.X, self.Y = X, Y
        theta = np.zeros(Y.shape[1] * X.shape[1])
        theta, j, ret = fmin_l_bfgs_b(self.cost_grad, theta, **self.fmin_args)
        if ret['warnflag'] != 0:
            warnings.warn('L-BFGS failed to converge')
        self.Theta = theta.reshape((Y.shape[1], X.shape[1]))

    def predict(self, X):
        P = np.exp(X.dot(self.Theta.T))
        P /= np.sum(P, axis=1)[:,None]

        return np.argmax(P, axis=1)


if __name__ == '__main__':
    def softmax():
        from data import load_svm

        X, y = load_svm('../data/iris.scale')
        Y = np.eye(3)[(y - 1 + 0.5).astype(np.int)]
        m = SoftmaxRegressionGD(lambda_=0)
        m.fit(X, Y)
        print np.mean(m.predict(X) == (y - 1 + 0.5).astype(np.int))

    softmax()

    def ex2():
        from data import load_txt

        X, y = load_txt('../data/ex2data1.txt')
        X = np.column_stack((np.ones(X.shape[0]), X))

        lr = LogisticRegressionGD(lambda_=0)
        lr.fit(X, y)
        print('Classifying [45, 85]:', lr.predict(np.array([1, 45, 85])))
        print('Training accuracy:', np.mean((lr.predict(X) > 0.5) == y))

        # plot data
        import matplotlib.pyplot as plt

        xs = [min(X[:, 1]), max(X[:, 1])]
        ys = [-(lr.theta[1] * x + lr.theta[0]) / lr.theta[2] for x in xs]
        plt.plot(xs, ys, '-')

        pos = y == 1
        neg = y == 0
        plt.plot(X[pos][:, 1], X[pos][:, 2], 'g+')
        plt.plot(X[neg][:, 1], X[neg][:, 2], 'ro')
        plt.show()

    def ex2_reg():
        from data import load_txt

        def map_feature(x, y):
            cols = []
            for sum in range(7):
                for i in range(sum + 1):
                    cols.append(x ** (sum - i) * y ** i)
            return np.column_stack(cols)

        X, y = load_txt('../data/ex2data2.txt')
        X_mapped = map_feature(X[:, 0], X[:, 1])

        lr = LogisticRegressionGD(lambda_=1)
        lr.fit(X_mapped, y)
        print('Training accuracy:', np.mean((lr.predict(X_mapped) > 0.5) == y))

        # plot data
        import matplotlib.pyplot as plt

        u = np.linspace(-1, 1.5, 50)
        v = np.linspace(-1, 1.5, 50)
        z = np.zeros((u.size, v.size))
        for i in range(u.size):
            for j in range(v.size):
                z[j, i] = map_feature(u[i], v[j]).dot(lr.theta)
        plt.contour(u, v, z, 0)

        pos = y == 1
        neg = y == 0
        plt.plot(X[pos][:, 0], X[pos][:, 1], 'g+')
        plt.plot(X[neg][:, 0], X[neg][:, 1], 'ro')
        plt.show()

    def ps1_1():
        import matplotlib.pyplot as plt

        y = np.loadtxt('../data/q1y.dat')
        X = np.loadtxt('../data/q1x.dat')
        X = np.column_stack((np.ones_like(y), X))

        lr = LogisticRegressionNR()
        lr.fit(X, y)
        print lr.theta

        pos = y == 0
        neg = y == 1

        plt.plot(X[pos,1], X[pos,2], 'rx')
        plt.plot(X[neg,1], X[neg,2], 'go')

        w = lr.theta
        x = np.array([X[:,1].min(), X[:,1].max()])
        y = (-w[1] * x - w[0]) / w[2]
        plt.plot(x, y)

        plt.show()
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.