# Hello world for linear regression in PyStan

This is meant to be a first notebook for exploring how to use PyStan to do:

* MAP estimation
* MCMC sampling using the No-U-Turn (NUTS) extension of Hamlitonian Monte Carlo

Meant to allow side-by-side comparison with alternative probabilistic programming packages like PyMC3.

Inspired by:

* Bob Carpenter's blog post:

https://statmodeling.stat.columbia.edu/2017/05/31/compare-stan-pymc3-edward-hello-world/

* The "hello world" example in the PyMC3 journal article:

Salvatier, Wiecki, and Fonnesbeck. "Probabilistic programming in Python using PyMC3". 
https://peerj.com/articles/cs-55.pdf


```python
import numpy as np
import pandas as pd
```


```python
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_style('whitegrid');
sns.set_context("notebook", font_scale=1.5, rc={"lines.linewidth": 2.0})

%matplotlib inline
```

# Import PyStan

If you need to install it first, using the `spr_2020s_env` conda environment just do:

```console
$ conda activate spr_2020s_env
$ conda install -c conda-forge pystan
```


```python
import pystan
```

# Model as MATH

Likelihood:
\begin{align}
 p(y_{1:N} | x_{1:N}, \alpha, \beta, \sigma) &= \text{NormalPDF}(\alpha + \beta_1 x_{n1} + \beta_2 x_{n2}, \sigma^2)
\end{align}

Priors on the variance of the likelihood:
\begin{align}
 p(\sigma) &= \text{HalfNormalPDF}(\sigma | 0, 1) = 
 \begin{cases}
     2 \cdot \text{NormalPDF}(\sigma | 0, 1) & \sigma > 0
\\
     0.0  & \sigma \leq 0
 \end{cases}
\end{align}


Priors on regression coefficients:
\begin{align}
 p(\alpha) &= \text{NormalPDF}(0, 10)
\\
 p(\beta_1) &= \text{NormalPDF}(0, 10)
\\
 p(\beta_2) &= \text{NormalPDF}(0, 10)
\\    
\end{align}

# Simulate data from the model


```python
# Initialize random number generator
np.random.seed(123)

# True regression coefficients
alpha = 1.0
beta = np.asarray([1.0, 1.337])

# True regression std dev
sigma = 0.3

# Size of dataset
N = 100

# Data Features (input variable)
x1_N = np.random.randn(N)
x2_N = np.random.randn(N)
x_N2 = np.hstack([x1_N[:,np.newaxis], x2_N[:,np.newaxis]])

# Simulate outcome variable
y_N = alpha + beta[0]*x1_N + beta[1]*x2_N + np.random.randn(N) * sigma

true_params = dict(
    alpha=alpha,
    beta=beta,
    sigma=sigma)
```

# Plot of observed data


```python
_, axgrid = plt.subplots(nrows=1, ncols=2, figsize=(8,3), sharex=True, sharey=True);

axgrid[0].plot(x1_N, y_N, 'k.'); axgrid[0].set_xlabel('x1'); axgrid[0].set_ylabel('y');
axgrid[1].plot(x2_N, y_N, 'b.'); axgrid[1].set_xlabel('x2'); axgrid[1].set_ylabel('y');
```

# Model specified in Stan


```python
model_code = '''
data {
  int N;
  vector[N] y;
  matrix[N, 2] x;
}
parameters {
  real alpha;
  vector[2] beta;
  real<lower=0> sigma; // constrain sigma > 0
}
model {
  alpha ~ normal(0, 10);
  beta ~ normal(0, 10);  // implicitly applied to each dim of vector beta
  sigma ~ normal(0, 1);  // implicitly takes into account the constraint: sigma > 0
  
  y ~ normal(alpha + x * beta, sigma);
}
'''
```


```python
# Create dictionary to specify concrete values
# for each parameter named in 'data' section of model

data_N10 = dict(N=10, y=y_N[:10], x=x_N2[:10])
```


```python
data_N100 = dict(N=100, y=y_N[:100], x=x_N2[:100])
```

# Setup: Compile the model


```python
# Compile the model
stan_model = pystan.StanModel(model_code=model_code)
```

    INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_9f0db3f7b821fb462af2f1e2557f82cf NOW.


# Setup: Prepare for plotting results


```python
def plot_traces(fit):
    fig, ax_grid = plt.subplots(nrows=2, ncols=2, sharex=True, sharey=True);
    for ax, var_name, dim in zip(
            ax_grid.flatten(), ['alpha', 'sigma', 'beta', 'beta'], [None, None, 0, 1]):
        var_samples_S = fit[var_name]
        true_val = true_params[var_name]

        if dim is not None:
            var_samples_S = var_samples_S[:, dim]
            true_val = true_val[dim]
            var_name = var_name + "[%d]" % dim

        S = var_samples_S.size
        ax.plot(var_samples_S, 'b.', markersize=4, alpha=0.3, label='samples')
        ax.plot(true_val * np.ones(S), 'r-', label='true_value')
        ax.set_ylabel(var_name)

    ax.legend(bbox_to_anchor=(1.1, 1.0))
    ax.set_ylim([0, 2.0]);
```

# Run MAP estimator for N=10 dataset


```python
map_N10 = stan_model.optimizing(
    data=data_N10,
    algorithm='LBFGS',
    init='0',
    verbose=True,
    iter=1000, seed=101)
```


```python
for key, arr in map_N10.items():
    print(key, arr)
```

    alpha 1.0693495791082472
    beta [1.00373535 1.44519637]
    sigma 0.2004478445079128



```python
for key, arr in true_params.items():
    print(key, arr)
```

    alpha 1.0
    beta [1.    1.337]
    sigma 0.3


# Run MAP estimator for N=100 dataset


```python
map_N100 = stan_model.optimizing(
    data=data_N100,
    algorithm='LBFGS',
    init='0',
    verbose=True,
    iter=1000, seed=101)
```


```python
for key, arr in map_N10.items():
    print(key, arr)
```

    alpha 1.0693495791082472
    beta [1.00373535 1.44519637]
    sigma 0.2004478445079128


# Run MCMC sampler for N=10 dataset


```python
fit_N10 = stan_model.sampling(
    data=data_N10,
    algorithm='NUTS',
    init='0',
    verbose=True,
    iter=1000, chains=1, warmup=500, thin=1, seed=101)
```


```python
plot_traces(fit_N10)
```

# Run MCMC sampler for N=100 dataset


```python
fit_N100 = stan_model.sampling(
    data=data_N100,
    algorithm='NUTS',
    init='0',
    verbose=True,
    iter=1000, chains=1, warmup=500, thin=1, seed=101)
```


```python
plot_traces(fit_N100)
```


```python

```

c312773d0cc31d202c066fc108724064c8cec92b
ipynb
Jupyter Notebook
notebooks/HelloWorldLinearRegression_Stan.ipynb
tufts-ml-courses/comp136-spr-20s-assignments-
c53cce8e376862eeef395aa0b55eca8b284a0115
MIT
2020-04-18T21:03:04.000Z
2020-04-18T21:03:04.000Z
notebooks/HelloWorldLinearRegression_Stan.ipynb
tufts-ml-courses/comp136-spr-20s-assignments-
c53cce8e376862eeef395aa0b55eca8b284a0115
MIT
notebooks/HelloWorldLinearRegression_Stan.ipynb
tufts-ml-courses/comp136-spr-20s-assignments-
c53cce8e376862eeef395aa0b55eca8b284a0115
MIT
2020-01-28T22:47:07.000Z
2020-04-12T23:56:18.000Z
Qwen/Qwen-72B
1. YES
2. YES