Pyglmnet

[//]: (----------------------------------------------------------------------)
class: center middle
# Pyglmnet
.center[
<img src="figures/pyglmnet.png" style="width: 700px;"/>
<div>
  <i class="fa fa-github"></i> github.com/glm-tools/pyglmnet
</div>
<div>
  <i class="fa fa-newspaper-o"></i> glm-tools.github.io/pyglmnet
</div>
]

.left[
<ul class="fa-ul">
  <li><i class="fa-li fa fa-github"></i>pavanramkumar</li>
  <li><i class="fa-li fa fa-twitter"></i>@desipoika</li>
  <li><i class="fa-li fa fa-home"></i>pavanramkumar.github.io</li>
</ul>
]
---

[//]: (----------------------------------------------------------------------)
# Why this package?
**Generalized linear models** encompass a broad class of models.
<img src="figures/glmpaper-00.png" style="width: 500px;"/>

**Elastic net** is an excellent algorithm for regularization.
<img src="figures/enetpaper-00.png" style="width: 500px;"/>

---

[//]: (----------------------------------------------------------------------)
# Why this package?

Fast implementation in widely popular **R** package **`glmnet`** for
**Generalized Linear Models**.

Scikit-learn has an implementation only for linear and logistic models.

|               | `normal`  | `logistic` | `poisson`  | `softplus` | `multinomial` | `cox`    |
|:--------------|:---------:|:----------:|:----------:|:----------:|:-------------:|:--------:|
|`R`            |     `x`   |  `x`       |  `x`       |  --        |  --           |  `x`     |
|`scikit-learn` |     `x`   |  `x`       |   --       |  --        |  --           |  --      |
|`pyglmnet`     |     `x`   |  `x`       |  `x`       |  `x`       |  `x`          |  --      |
---

[//]: (----------------------------------------------------------------------)
# @KordingLab, Northwestern

.center[
<img src="figures/pidata.png" style="width: 800px;"/>
]

---
[//]: (----------------------------------------------------------------------)
# Machine Learning for Brain and Behavior
.center[
<img src="figures/about-00.png" style="width: 600px;"/>
]
---

# We often want to model spike counts
.center[
<img src="figures/about-01.png" style="width: 600px;"/>
]
---

# Count statistics abound in the real world
<img src="figures/githubstars-01.png" style="width: 700px;"/>

---

[//]: (----------------------------------------------------------------------)
# Example neurons

---

# Linear vs. Poisson Models
<img src="figures/neurons-03.png" style="width: 600px;"/>

--
<img src="figures/neurons-04.png" style="width: 600px;"/>

See [Linear vs. Poisson notebook](https://nbviewer.jupyter.org/github/pavanramkumar/pydata-chicago-2016/blob/master/notebooks/Linear_vs_Poisson.ipynb).
---

[//]: (----------------------------------------------------------------------)
# Linear Regression

### Predict
- a target variable, `$y$`
- as a linear combination of `$p$` predictors, `$\mathbf{x} = [x_{1}, \dots, x_{p}]$`.
- from a set of `$n$` observations `$\left\{\mathbf{x}_i, y_i\right\}$`.

$$
\hat{y} = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p
$$

- In matrix form:

$$
\mathbf{\hat{y}} = \beta_0 + \mathbf{X} \mathbf{\beta}
$$

where `$\mathbf{\hat{y}} = [\hat{y}_1, \dots, \hat{y}_n]$`
and `$\mathbf{X} = [\mathbf{x}_1^T, \dots, \mathbf{x}_n^T]$`.

```python
def predict(X, beta0, beta):
    y = beta0 + np.dot(X, beta)
    return y
```
---

[//]: (----------------------------------------------------------------------)
# Linear Regression

### Mean Squared Loss
- To fit the parameters `$\left\{\beta_0, \beta\right\}$`,
we minimize the mean squared loss:
$$
J = \frac{1}{2n} \sum_i (y_i - \hat{y}_i)^2
= \frac{1}{2n} \sum_i (y_i - (\beta_0 + \mathbf{x}_i^T \mathbf{\beta}))^2
$$

- In matrix form:
$$
J = \frac{1}{2n} (\mathbf{y} - (\beta_0 + \mathbf{X} \mathbf{\beta}))^T
    (\mathbf{y} - (\beta_0 + \mathbf{X} \mathbf{\beta}))
$$

```python
def loss(X, y, beta0, beta):
    n = y.shape[0]
    yhat = predict(X, beta0, beta)
    err = y - yhat
    J = 1. / (2. * n) * np.dot(err.T , err)
    return J
```
---

[//]: (----------------------------------------------------------------------)
# Linear Regression

### Negative Log Likelihood Loss

- Real data is noisy. Let's assume random Gaussian noise.
$$
\mathbf{y} = \beta_0 + \mathbf{X}\beta + \epsilon; \epsilon \sim \mathcal{N}(\mathbf{0}, \sigma^2\mathbf{I})
$$

- The likelihood of observing the target variables, given the predictors is:
$$
P(\mathbf{y} | \mathbf{X}; \beta_0, \beta) \propto \prod_i \exp\left\[-\frac{1}{2} (y_i - \beta_0 - \mathbf{x}_i^T \beta)^2\right\]
$$

- The negative log-likelihood function is identical to the mean squared loss!
$$
-\mathcal{L}(\beta_0, \beta) = \frac{1}{2n}\sum_i (y_i - \hat{y}_i)^2
$$

???
Thus, we can write out the loss as negative log-likelihood for any distribution!
---

[//]: (----------------------------------------------------------------------)
# Generalized Linear Models

Recall that for linear regression
$$
\mathbf{y} = \beta_0 + \mathbf{X}\beta + \epsilon; \epsilon \sim \mathcal{N}(\mathbf{0}, \sigma^2\mathbf{I})
$$

Thus, `$\mathbf{y}$` is a multivariate normal distribution with mean function `$\mu$`.
$$
\mathbf{y} \sim \mathcal{N}(\mathbf{\mu}, \sigma^2\mathbf{I}); \mathbf{\mu} = \beta_0 + \mathbf{X}\beta
$$

This is a special case of GLMs:
$$
\mathbf{y} \sim \text{ExpFamily}(\mathbf{\mu}, \Theta); \mathbf{\mu} = q(\beta_0 + \mathbf{X}\beta)
$$

where
- `$q(.)$` is a pointwise nonlinearity, and
- `$\text{ExpFamily}(\mathbf{\mu}, \Theta)$` is any pdf in the exponential family,
with mean function `$\mathbf{\mu}$` and parameters `$\Theta$`.

???
Using this paramterization, GLMs allow us to model target variables `$\mathbf{y}$`
that take on different ranges (e.g. binary variables, counts)
and distributions, using real-valued predictors `$\mathbf{X}$`.
---

[//]: (----------------------------------------------------------------------)
# Example: Poisson Regression

Consider the prediction problem where `$y_i$`'s are count variables.
`$\mathbf{y} = \{y_i\}$` where `$y_i \in \{0, 1, 2, \dots \}$`.

#### Assumption 1: Nonlinearity
We set `$q(z) = \exp(z)$`,
with domain `$(-\infty, \infty)$`, and range `$[0, \infty)$`.

```python
def predict(X, beta0, beta):
    y = np.exp(beta0 + np.dot(X, beta))
    return y
```

#### Assumption 2: Noise distribution
We assume that `$y_i$`'s are Poisson random variables with mean `$\lambda$`.
We model this as
`$\lambda \sim \hat{y}_i = q(\beta_0 + \mathbf{x}_i^T\beta)$`.

---

[//]: (----------------------------------------------------------------------)
# Example: Poisson Regression

Let's write out the negative log-likelihood of the data given these assumptions.
$$
P(y_i = k) = \frac{e^{-\lambda} \lambda^k}{k!}
$$

$$
P(\mathbf{y} | \mathbf{X}; \beta_0, \beta) \propto
\prod_i e^{-\lambda_i} (\lambda_i)^{y_i}
$$

$$
-\mathcal{L}(\beta_0, \beta) = -\Bigg[\sum_i y_i log(\hat{y}_i) -\hat{y}_i \Bigg]
$$

This is known as the Poisson loss.
```python
def loss(X, y, beta0, beta):
    yhat = predict(X, beta0, beta)
    J = -np.sum(y * np.log(yhat) - yhat)
    return J
```
---

[//]: (----------------------------------------------------------------------)
# Example: Logistic Regression

Consider the binary classification scenario, where `$y_i$`'s are binary random variables.
`$\mathbf{y} = \{y_i\}$` where `$y_i \in \{0, 1\}$`.

#### Assumption 1: Nonlinearity
We set `$q(z) = \sigma(z) = \frac{1}{1 + e^{-z}}$`,
which is the sigmoid function, with domain `$(-\infty, \infty)$`, and range `$[0, 1]$`.

```python
def sigmoid(z):
    s = np.exp(z) / (1 + np.exp(z))
    return s

def predict(X, beta0, beta):
    y = sigmoid(beta0 + np.dot(X, beta))
    return y
```
---

[//]: (----------------------------------------------------------------------)
# Example: Logistic Regression

#### Assumption 2: Noise distribution
We assume that `$y_i$`'s are Bernoulli random variables with mean `$p$`.
`$P(y_i = 1) = p; P(y_i = 0) = 1-p$`.
We model `$p \sim \hat{y}_i = q(\beta_0 + \mathbf{x}_i^T\beta)$`.

Let's write out the negative log-likelihood of the data given these assumptions.

$$
P(\mathbf{y} | \mathbf{X}; \beta_0, \beta) =
\prod_i p^{y_i} (1-p)^{(1-y_i)}
$$

$$
-\mathcal{L}(\beta_0, \beta) = -\Bigg[\sum_i y_i log(\hat{y}_i) + (1-y_i) log(1-\hat{y}_i) \Bigg]
$$

This is known as the cross entropy loss.

```python
def loss(X, y, beta0, beta):
    yhat = predict(X, beta0, beta)
    J = -np.sum(y * np.log(yhat) + \
               (1 - y) * np.log(1 - yhat))
    return J
```
---

[//]: (----------------------------------------------------------------------)
# Optimization

Recalling our calculus 101, we minimize the loss `$J$` by taking the derivative
(aka gradient) and setting it to zero.

#### Linear Regression
.left[`$
J = \frac{1}{2n} (\mathbf{y} - \mathbf{\hat{y}})^T
    (\mathbf{y} - \mathbf{\hat{y}})
$`]

.left[`$
\frac{\partial J}{\partial \beta_0}
=  -\frac{1}{n}(\mathbf{y} - \mathbf{\hat{y}})
= -\frac{1}{n}(\mathbf{y} - \beta_0 - \mathbf{X}\beta)
= 0
$`]
.left[`$
\frac{\partial J}{\partial \beta}
=  -\frac{1}{n}\mathbf{X}^T(\mathbf{y} - \mathbf{\hat{y}})
= -\frac{1}{n}\mathbf{X}^T(\mathbf{y} - \beta_0 - \mathbf{X}\beta)
= 0
$`]

.left[
```python
def grad_loss(X, y, beta0, beta):
    n = y.shape[0]
    yhat = predict(X, beta0, beta)
    err = y - yhat
    grad_beta0 = -1. / n * err
    grad_beta = -1. / n * np.dot(X.T, err)
    return grad_beta0, grad_beta
```
]
---

[//]: (----------------------------------------------------------------------)
# Optimization

Recalling our calculus 101, we minimize the loss `$J$` by taking the derivative
(aka gradient) and setting it to zero.

#### Logistic Regression
.left[`$
J = -\frac{1}{n} \left\{ \mathbf{y}^T \log\left\{\mathbf{\hat{y}}\right\} +
          (\mathbf{1} - \mathbf{y})^T \log\left\{\mathbf{1} - \mathbf{\hat{y}}\right\}
          \right\}
$`]

.left[`$
\frac{\partial J}{\partial \beta_0}
=  -\frac{1}{n}(\mathbf{y} - \mathbf{\hat{y}})
= -\frac{1}{n}(\mathbf{y} - \mathbf{1}./ (\mathbf{1} + \exp(-\beta_0 - \mathbf{X}\beta)))
= 0
$`]
.left[`$
\frac{\partial J}{\partial \beta}
=  -\frac{1}{n}\mathbf{X}^T(\mathbf{y} - \mathbf{\hat{y}})
= -\frac{1}{n}\mathbf{X}^T(\mathbf{y} - \mathbf{1}./ (\mathbf{1} + \exp(-\beta_0 - \mathbf{X}\beta)))
= 0
$`]

[//]: (----------------------------------------------------------------------)
# Gradient Descent

We cannot always solve this system of equations in closed form.
But we can arrive at a unique solution by guessing a starting point and
iteratively taking small steps against the gradient.
This is called gradient descent.