Jure Žbontar avatar Jure Žbontar committed 7a2b7e7

Neural networks for regression.

Comments (0)

Files changed (2)

mlclass/linear_regression.py

         theta = np.zeros(X.shape[1])
         self.theta, j, ret = fmin_l_bfgs_b(self.cost_grad, theta, args=(X, y),
             **self.fmin_args)
-        assert ret['warnflag'] == 0
 
     def predict(self, X):
         return X.dot(self.theta)

mlclass/neural_networks.py

 import numpy as np
-import scipy.sparse
+import scipy.sparse as sp
 from scipy.optimize import fmin_l_bfgs_b
 
 class NeuralNetwork:
         eps = np.sqrt(6) / np.sqrt(i0 + i2)
         initial_params = np.random.randn(n_params) * 2 * eps - eps
 
-        self.X = self.append_ones(X)
+        self.X = append_ones(X)
         self.y = y
 
         params, _, _ = fmin_l_bfgs_b(self.cost_grad, initial_params, **self.fmin_args)
     def predict(self, X):
         m, n = X.shape
         
-        a2 = sigmoid(self.append_ones(X).dot(self.Theta1.T))
+        a2 = sigmoid(append_ones(X).dot(self.Theta1.T))
         a3 = sigmoid(np.column_stack((np.ones(m), a2)).dot(self.Theta2.T))
 
         return a3
 
-    def append_ones(self, 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 append_ones(X):
+    if sp.issparse(X):
+        return sp.hstack((np.ones((X.shape[0], 1)), X)).tocsr()
+    else:
+        return np.hstack((np.ones((X.shape[0], 1)), X))
 
 def sigmoid(x):
     return 1.0 / (1.0 + np.exp(-x))
     sx = sigmoid(x)
     return sx * (1 - sx)
 
+class MLPRegressor:
+    def __init__(self, layers, lambda_=1, callback=None, **fmin_args):
+        self.layers = layers
+        self.lambda_ = lambda_
+        self.callback = callback
+        self.fmin_args = fmin_args
+
+    def unfold_params(self, params):
+        i0, i1, i2 = self.layers
+
+        i = (i0 + 1) * i1
+
+        Theta1 = params[:i].reshape((i1, i0 + 1))
+        Theta2 = params[i:].reshape((i2, i1 + 1))
+
+        return Theta1, Theta2
+
+    def cost_grad(self, params, X, y):
+        Theta1, Theta2 = self.unfold_params(params)
+
+        if self.callback:
+            self.Theta1 = Theta1
+            self.Theta2 = Theta2
+            self.callback(self)
+
+        # Feedforward Propagation
+        m, n = X.shape
+
+        a1 = append_ones(X)
+        z2 = a1.dot(Theta1.T)
+        a2 = append_ones(sigmoid(z2))
+        z3 = a2.dot(Theta2.T)
+        a3 = z3
+
+        # Cost
+        J = np.sum((a3 - y)**2)
+
+        t1 = Theta1.copy()
+        t1[:, 0] = 0
+        t2 = Theta2.copy()
+        t2[:, 0] = 0
+
+        # regularization
+        reg = np.dot(t1.flat, t1.flat)
+        reg += np.dot(t2.flat, t2.flat)
+        J += float(self.lambda_) * reg / 2.0
+
+        J /= m
+
+        # Grad
+        d3 = 2 * (a3 - y)
+        d2 = d3.dot(Theta2)[:, 1:] * sigmoid_gradient(z2)
+
+        D2 = a2.T.dot(d3).T / m
+        D1 = a1.T.dot(d2).T / m
+
+        # regularization
+        D2 += t2 * (float(self.lambda_) / m)
+        D1 += t1 * (float(self.lambda_) / m)
+
+        return J, np.hstack((D1.flat, D2.flat))
+
+    def fit(self, X, y):
+        i0, i1, i2 = self.layers
+
+        m, n = X.shape
+        n_params = i1 * (i0 + 1) + i2 * (i1 + 1)
+        eps = np.sqrt(6) / np.sqrt(i0 + i2)
+        initial_params = np.random.randn(n_params) * 2 * eps - eps
+
+        params, _, _ = fmin_l_bfgs_b(self.cost_grad, initial_params, args=(X, y),
+            **self.fmin_args)
+        self.Theta1, self.Theta2 = self.unfold_params(params)
+
+    def predict(self, X):
+        m, n = X.shape
+        
+        a2 = sigmoid(append_ones(X).dot(self.Theta1.T))
+        a3 = append_ones(a2).dot(self.Theta2.T)
+
+        return a3
+
 if __name__ == '__main__':
+    def mlp_regressor():
+        from sklearn.datasets import load_boston
+
+        data = load_boston()
+        X = data['data']
+        y = data['target']
+
+        layers = [X.shape[1], 20, 1]
+        m = MLPRegressor(layers, 1)
+        m.fit(X, y[:,None])
+        p1 = m.predict(X).ravel()
+
+        from mlclass.linear_regression import LinearRegressionGD
+
+        m = LinearRegressionGD(1)
+        m.fit(X, y)
+        p2 = m.predict(X)
+
+        print np.column_stack((p1, p2, y))
+        
+
+    mlp_regressor()
+
     def ex3():
         from scipy.io import loadmat
 
 
         nn = NeuralNetwork([400, 25, 10], iprint=1, maxfun=10)
         nn.fit(X, Y)
-
-    test()
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.