# Machine Learning Implementation

## Imports


```python
import json

import numpy as np
import pandas as pd
import plotly.offline as py
from plotly import graph_objects as go
```

## Logistic regression

### The maths

The logistic model aims to predict the discrete y variable a.k.a the target variable (e.g. whether something will happen) based on a collection of features. It does this by transforming a linear combination of the features into a curve and fitting this curve to the data.

The curve used in logistic regression is the sigmoid function

$$
\sigma(x) = \frac{1}{1+e^{-x}}
$$

Define y as

$$
\begin{align}
\hat{y} &= h_{\boldsymbol{\beta}}(\mathbf{x})\\
\hat{y}&= \sigma\left(\beta_0x_0+\cdots+\beta_nx_n\right)\quad &n\in \mathbb{N},x_0=1 \\
\hat{y}&=\sigma\left(\sum^{n}_{i=0}\beta_ix_i\right) \\
\hat{y}&=\sigma\left(\mathbf{\boldsymbol{\beta}^Tx}\right)\quad&\boldsymbol{\beta},\mathbf{x}\in\mathbb{R}^{n\times1}\\
\hat{y}&=\sigma\left(\boldsymbol{\beta}^T\mathbf{x}\right)
\end{align}
$$

notice

$$
\hat{y} = \frac{1}{1+e^{-\boldsymbol{\beta}^T\mathbf{x}}}
$$

so

$$
\begin{align}
\hat{y} + \hat{y}e^{-\boldsymbol{\beta}^T\mathbf{x}} &= 1\\
\hat{y}e^{-\boldsymbol{\beta}^T\mathbf{x}} &= 1 - \hat{y}\\
\frac{\hat{y}}{1 - \hat{y}} &= e^{\boldsymbol{\beta}^T\mathbf{x}}\\
\ln\left(\frac{\hat{y}}{1 - \hat{y}}\right)&=\boldsymbol{\beta}^T\mathbf{x}
\end{align}
$$

This above is the logit form of logistic regression. We model the logit as a linear combination of the x variables

We define the cost function as follows for each y and corresponding x

$$
\begin{align}
J(\mathbf{x})
&= \begin{cases}
-\log\left(h_{\boldsymbol{\beta}}(\mathbf{x})\right) &\text{if y=1}\\
-\log\left(1-h_{\boldsymbol{\beta}}(\mathbf{x})\right) &\text{if y=0}\\
\end{cases}
\end{align}
$$

$$
\begin{align}
J(\mathbf{x})
&= -\frac{1}{m}\sum_{j=1}^my^j\log\left(h_{\boldsymbol{\beta}}(\mathbf{x}^j)\right)
+(1-y^j)\log\left(1-h_{\boldsymbol{\beta}}(\mathbf{x}^j)\right)\\
&= -\frac{1}{m}\sum_{j=1}^my^j\log\left(\frac{1}{1+e^{-\boldsymbol{\beta}^T\mathbf{x}}}\right)
+(1-y^j)\log\left(1-\frac{1}{1+e^{-\boldsymbol{\beta}^T\mathbf{x}}}\right)\\
&= -\frac{1}{m}\sum_{j=1}^my^j\log\left(\frac{1}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}\right)
+(1-y^j)\log\left(1-\frac{1}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}\right)
\end{align}
$$

note

$$
\begin{align}
h_{\boldsymbol{\beta}}(\mathbf{x}^j)&=\frac{1}{1+e^{-\boldsymbol{\beta}^T\mathbf{x}^j}}\\
&=\frac{1}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}
\end{align}
$$

so

$$
\begin{align}
\frac{\partial h}{\partial \beta_k} &= -\left(1+e^{-\sum^{n}_{i=0}\beta_ix_i}\right)^{-2}e^{-\sum^{n}_{i=0}\beta_ix_i} (-x_k^j)\\
&=\frac{1}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}
\frac{-e^{-\sum^{n}_{i=0}\beta_ix_i} (-x_k^j)}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}\\
&=\frac{1}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}
\frac{(1-1-e^{-\sum^{n}_{i=0}\beta_ix_i})(-x_k^j)}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}\\
&=\frac{1}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}
\left(
\frac{1}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}-
\frac{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}
\right)(-x_k^j)\\
&=\frac{1}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}
\left(
\frac{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}-
\frac{1}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}
\right)(x_k^j)\\
&=\frac{1}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}
\left(
1-
\frac{1}{1+e^{-\sum^{n}_{i=0}\beta_ix_i}}
\right)(x_k^j)\\
&=h_{\boldsymbol{\beta}}(\mathbf{x}^j)(1-h_{\boldsymbol{\beta}}(\mathbf{x}^j))x_k^j
\end{align}
$$

We need to differentiate the cost function i.e. find the gradient

$$
\begin{align}
\frac{\partial J}{\partial\beta_k}\left(\boldsymbol{\beta}\right) 
&=\frac{\partial}{\partial\beta_k}\left(
-\frac{1}{m}\sum_{j=1}^my^j\log\left(h_{\boldsymbol{\beta}}(\mathbf{x}^j)\right)
+(1-y^j)\log\left(1-h_{\boldsymbol{\beta}}(\mathbf{x}^j)\right)
\right)\\
&=-\frac{1}{m}\sum_{j=1}^m\frac{y^j}{h_{\boldsymbol{\beta}}(\mathbf{x}^j)}\frac{\partial h}{\partial \beta_k}
+\frac{-(1-y^j)}{1-h_{\boldsymbol{\beta}}(\mathbf{x}^j)}\frac{\partial h}{\partial \beta_k}\\
&=-\frac{1}{m}\sum_{j=1}^m\frac{y^j}{h_{\boldsymbol{\beta}}(\mathbf{x}^j)}
h_{\boldsymbol{\beta}}(\mathbf{x}^j)(1-h_{\boldsymbol{\beta}}(\mathbf{x}^j))x_k^j
+\frac{-(1-y^j)}{1-h_{\boldsymbol{\beta}}(\mathbf{x}^j)}
h_{\boldsymbol{\beta}}(\mathbf{x}^j)(1-h_{\boldsymbol{\beta}}(\mathbf{x}^j))x_k^j\\
&=-\frac{1}{m}\sum_{j=1}^my^j(1-h_{\boldsymbol{\beta}}(\mathbf{x}^j))x_k^j
-(1-y^j)
h_{\boldsymbol{\beta}}(\mathbf{x}^j)x_k^j\\
&=\frac{1}{m}\sum_{j=1}^m
\left(h_{\boldsymbol{\beta}}(\mathbf{x}^j)-y^j\right)x_k^j
\end{align}
$$

hence

$$
\nabla_{\boldsymbol{\beta}} J
=
\begin{bmatrix}
       \frac{\partial J}{\partial\beta_1} \\
       \vdots \\
       \frac{\partial J}{\partial\beta_n}
\end{bmatrix}
=
\begin{bmatrix}
       \frac{1}{m}\sum_{j=1}^m
            \left(h_{\boldsymbol{\beta}}(\mathbf{x}^j)-y^j\right)x_1^j\\
       \vdots \\
       \frac{1}{m}\sum_{j=1}^m
           \left(h_{\boldsymbol{\beta}}(\mathbf{x}^j)-y^j\right)x_n^j
\end{bmatrix}
$$

Define the design matrix and column representation of y. Here each row of X and y are training examples hence there are m rows

$$\mathbf{X}\in\mathbb{R}^{m\times n},
\quad \mathbf{y}\in\mathbb{R}^{m\times 1}
$$

$$
\mathbf{X}=\begin{bmatrix}
       \dots & (\mathbf{x}^1)^T & \dots\\
       \dots & (\mathbf{x}^2)^T & \dots\\
       \dots & \vdots  & \dots\\
       \dots & (\mathbf{x}^m)^T & \dots
\end{bmatrix}\quad
\mathbf{y}=\begin{bmatrix}
    y_1\\y_2\\\vdots\\y_m
\end{bmatrix}
$$

$$
\begin{align}
\nabla_{\boldsymbol{\beta}} J
=
\begin{bmatrix}
       \frac{1}{m}\sum_{j=1}^m
            \left(h_{\boldsymbol{\beta}}(\mathbf{x}^j)-y^j\right)x_1^j\\
       \vdots \\
       \frac{1}{m}\sum_{j=1}^m
           \left(h_{\boldsymbol{\beta}}(\mathbf{x}^j)-y^j\right)x_n^j
\end{bmatrix}
=
\frac{1}{m}
\begin{bmatrix}
       \sum^{n}_{i=0}h_{\boldsymbol{\beta}}(\mathbf{x}^j)^jx^j_1\\
       \vdots \\
       \sum^{n}_{i=0}h_{\boldsymbol{\beta}}(\mathbf{x}^j)x^j_n
\end{bmatrix}
-
\frac{1}{m}
\begin{bmatrix}
       \sum^{m}_{j=1}y^jx^j_1\\
       \vdots \\
       \sum^{m}_{j=1}y^jx^j_n\\
\end{bmatrix}
\end{align}
$$

$$
h_{\boldsymbol{\beta}}(\mathbf{x}^j) = \sigma({\mathbf{x}^j}^T\boldsymbol{\beta})
$$

so

$$
\begin{align}
\nabla_{\boldsymbol{\beta}} J
&=\frac{1}{m}\left(
\mathbf{X}^T\sigma(\mathbf{X}\mathbf{\boldsymbol{\beta}})-\mathbf{X}^T\mathbf{y}
\right)\\
&=\frac{1}{m}\mathbf{X}^T\left(
\sigma(\mathbf{X}\mathbf{\boldsymbol{\beta}})-\mathbf{y}
\right)\\
&=\frac{1}{m}\mathbf{X}^T\left(
\mathbf{\hat{y}}-\mathbf{y}
\right)
\end{align}
$$

where

$$
\mathbf{\hat{y}} = \sigma(\mathbf{X}\mathbf{\boldsymbol{\beta}})
$$

We could have derived the same thing using matrix calculus

### Example sigmoid

The curve used in logistic regression is the sigmoid function

$$
\sigma(x) = \frac{1}{1+e^{-x}}
$$


```python
sigmoid_fig = go.FigureWidget()
demo_x = np.arange(-10,10,0.1)
demo_y = 1 / (1 + np.exp(-demo_x))
sigmoid_fig.add_scatter(
    x=demo_x,
    y=demo_y)
sigmoid_fig.layout.title = 'Sigmoid Function'
sigmoid_fig
```


    FigureWidget({
        'data': [{'type': 'scatter',
                  'uid': 'fc81bfbe-f9f8-419c-9923-4d127962d5e2',
     …


### Make fake data


```python
m = 100
x0 = np.ones(shape=(m, 1))
x1 = np.linspace(0, 10, m).reshape(-1, 1)
X = np.column_stack((x0, x1))

# let y = 0.5 * x + 1 + epsilon
epsilon =  np.random.normal(scale=2, size=(m, 1))
y = x1 + epsilon
y = (y > 5).astype(int)
```


```python
fig = go.FigureWidget()
fig = fig.add_scatter(
    x=X[:,1],
    y=y[:,0],
    mode='markers',
    name='linear data + noise')
fig
```


    FigureWidget({
        'data': [{'mode': 'markers',
                  'name': 'linear data + noise',
                  'ty…


### Logistic regression class


```python
import json

import numpy as np


class LogisticRegression():

    def __init__(self, learning_rate=0.05):
        """  
        Logistic regression model

        Parameters:
        ----------
        learning_rate: float, optional, default 0.05
            The learning rate parameter controlling the gradient descent
            step size
        """
        self.learning_rate = learning_rate
        print('Creating logistic model instance')

    def __repr__(self):
        return (
            f'<LogisticRegression '
            f'learning_rate={self.learning_rate}>')

    def fit(self, X, y, n_iter=1000):
        """  
        Fit the logistic regression model

        Updates the weights with n_iter iterations of batch gradient
        descent updates

        Parameters:
        ----------
        X: numpy.ndarray
            Training data, shape (m samples, (n - 1) features + 1)
            Note the first column of X is expected to be ones (to allow 
            for the bias to be included in beta)
        y: numpy.ndarray
            Target values - class label {0, 1}, shape (m samples, 1)
        n_iter: int, optional, default 1000
            Number of batch gradient descent steps
        """
        m, n = X.shape
        print(f'fitting with m={m} samples with n={n-1} features\n')
        self.beta = np.zeros(shape=(n, 1))
        self.costs = []
        self.betas = [self.beta]
        for iteration in range(n_iter):
            y_pred = self.predict_proba(X)
            cost = (-1 / m) * (
                (y.T @ np.log(y_pred)) +
                ((np.ones(shape=y.shape) - y).T @ np.log(
                    np.ones(shape=y_pred.shape) - y_pred))
            )
            self.costs.append(cost[0][0])
            gradient = (1 / m) * X.T @ (y_pred - y)
            self.beta = self.beta - (
                self.learning_rate * gradient)
            self.betas.append(self.beta)

    def predict_proba(self, X):
        """  
        Predicted probability values for class 1

        Note this is calculated as the sigmoid of the linear combination
        of the feature values and the weights.

        Parameters:
        ----------
        X: numpy.ndarray
            Training data, shape (m samples, (n - 1) features + 1)
            Note the first column of X is expected to be ones (to allow 
            for the bias to be included in beta)

        Returns:
        -------
        numpy.ndarray:
            Predicted probability of samples being in class 1
        """        
        y_pred = self.sigmoid(X @ self.beta)
        return y_pred

    def predict(self, X, descision_prob=0.5):
        """  
        Predict the class values from sample X feature values

        Parameters:
        ----------
        X: numpy.ndarray
            Training data, shape (m samples, (n - 1) features + 1)
            Note the first column of X is expected to be ones (to allow 
            for the bias to be included in beta)

        Returns:
        -------
        numpy.ndarray:
            Prediceted class values, shape (m samples, 1)
        """
        y_pred = self.sigmoid(X @ self.beta)
        return (y_pred > descision_prob) * 1

    def sigmoid(self, x):
        """  
        Sigmoid function

        f(x) = 1 / (1 + e^(-x))

        Parameters:
        ----------
        x: numpy.ndarray

        Returns:
        -------
        numpy.ndarray:
            sigmoid of x, values in (0, 1)
        """        
        return 1 / (1 + np.exp(-x))

```


```python
logistic_regression = LogisticRegression()
logistic_regression.fit(X, y)
```

    Creating logistic model instance
    fitting with m=100 samples with n=1 features
    



```python
example_X = np.array([[1,1],[1,4],[1,7]])
logistic_regression.predict(example_X)
```




    array([[0],
           [0],
           [1]])



### Plot the best fit


```python
fig = fig.add_scatter(
    x=X[:,1], 
    y=logistic_regression.predict_proba(X)[:,0],
    mode='markers',
    name='logistic best fit')
fig
```


    FigureWidget({
        'data': [{'mode': 'markers',
                  'name': 'linear data + noise',
                  'ty…


### Plot the cost function


```python
# Haven't got round to this yet - see linear regression for an example error 
# surface decent.
```

## Logisitc regression - Titanic example

### Load data


```python
X_train = pd.read_feather('../data/titanic/processed/X_train.feather')
X_test = pd.read_feather('../data/titanic/processed/X_test.feather')
y_train = pd.read_feather('../data/titanic/processed/y_train.feather')
y_test = pd.read_feather('../data/titanic/processed/y_test.feather')
```

### Train model


```python
titanic_logistic_model = LogisticRegression()
```

    Creating logistic model instance



```python
titanic_logistic_model.fit(X=X_train.values, y=y_train.values, n_iter=4000)
```

    fitting with m=712 samples with n=29 features
    


### Plot the cost


```python
titanic_cost_fig = go.FigureWidget()

titanic_cost_fig.add_scatter(
    x=list(range(len(titanic_logistic_model.costs))),
    y=titanic_logistic_model.costs,
)

titanic_cost_fig.layout.title = 'Cost Vs gradient descent iterations'
titanic_cost_fig.layout.xaxis.title = 'Iterations'
titanic_cost_fig.layout.yaxis.title = 'Cost'
titanic_cost_fig
```


    FigureWidget({
        'data': [{'type': 'scatter',
                  'uid': 'a3a5dfba-fddd-428e-87a9-4a7f4461bffc',
     …


### Error analysis


```python
y_pred = titanic_logistic_model.predict(X_test.values)
test_accuracy = (y_pred == y_test.values).sum() / len(y_pred)

print(f'Test accuracy is with my implementation {test_accuracy:.2%}')
```

    Test accuracy is with my implementation 79.89%


## End

9a744cf238de12ba8277625f6c33ad4e05fe912d
ipynb
Jupyter Notebook
notebooks/logistic_regression.ipynb
simonwardjones/machine_learning
1e92865bfe152acaf0df2df8f11a5f51833389a9
MIT
2020-05-13T13:21:14.000Z
2021-03-29T08:11:22.000Z
notebooks/logistic_regression.ipynb
simonwardjones/machine_learning
1e92865bfe152acaf0df2df8f11a5f51833389a9
MIT
notebooks/logistic_regression.ipynb
simonwardjones/machine_learning
1e92865bfe152acaf0df2df8f11a5f51833389a9
MIT
2020-06-20T12:50:23.000Z
2020-08-19T13:29:01.000Z
Qwen/Qwen-72B
1. YES
2. YES