# Introduction to probabilistic programming with PyMC3


```python
import os,sys
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import scipy.stats as st

%matplotlib inline
%precision 4
plt.style.use('bmh')
```

## Coin toss example

We are looking to determine if our coin is biased.  We can investigate the fairness of a coin in several ways.  One useful distribution is the Binomial.  It is parameterized by the number of trials ($n$) and the success probability in each trial ($p$).


```python
from IPython.display import Image
Image(filename='../binomial.png')
```


```python
## lets start with some data
n = 100
pcoin = 0.62 # actual value of p for coin
results = st.bernoulli(pcoin).rvs(n)
h = sum(results)
print("We observed %s heads out of %s"%(h,n))

## the expected distribution
p = 0.5
rv = st.binom(n,p)
mu = rv.mean()
sd = rv.std()
print("The expected distribution for a fair coin is mu=%s, sd=%s"%(mu,sd))

## p-value by simulation
nsamples = 100000
xs = np.random.binomial(n, p, nsamples)
print("Simulation p-value - %s"%(2*np.sum(xs >= h)/(xs.size + 0.0)))

## p-value by binomial test
print("Binomial test - %s"%st.binom_test(h, n, p))

## MLE
print("Maximum likelihood %s"%(np.sum(results)/float(len(results))))

## bootstrap
bs_samples = np.random.choice(results, (nsamples, len(results)), replace=True)
bs_ps = np.mean(bs_samples, axis=1)
bs_ps.sort()
print("Bootstrap CI: (%.4f, %.4f)" % (bs_ps[int(0.025*nsamples)], bs_ps[int(0.975*nsamples)]))
```

    We observed 59 heads out of 100
    The expected distribution for a fair coin is mu=50.0, sd=5.0
    Simulation p-value - 0.08668
    Binomial test - 0.0886260801141
    Maximum likelihood 0.59
    Bootstrap CI: (0.4900, 0.6900)


### Lets perform inference with PyMC3


```python
print("n = %s"%n)
print("h = %s"%h)
alpha = 2
beta = 2

niter = 1000
with pm.Model() as model: # context management
    # define priors
    p = pm.Beta('p', alpha=alpha, beta=beta)

    # define likelihood
    y = pm.Binomial('y', n=n, p=p, observed=h)

    # inference
    start = pm.find_MAP() # Use MAP estimate (optimization) as the initial state for MCMC
    step = pm.Metropolis() # Have a choice of samplers
    trace = pm.sample(niter, step, start, random_seed=123, progressbar=True)
```

    n = 100
    h = 59


    Applied logodds-transform to p and added transformed p_logodds_ to model.
    100%|██████████| 1000/1000 [00:00<00:00, 8636.83it/s]



```python
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)

ax.hist(trace['p'], 15, histtype='step', normed=True, label='post');
x = np.linspace(0, 1, 100)
ax.plot(x, st.beta.pdf(x, alpha, beta), label='prior');
ax.legend(loc='best');
```

## Estimating the mean and standard deviation of a normal

One of the amazing things about probabilistic programming in general is that to evaluate a complex model it is just an extension of what we would do with simple models---we go about it using the same process.

$$
X \sim \mathcal{N}(\mu,\sigma^{2})
$$


```python
# generate observed data
N = 100
_mu = np.array([10])
_sigma = np.array([2])
y = np.random.normal(_mu, _sigma, N)

niter = 1000
with pm.Model() as model:
    # define priors
    mu = pm.Uniform('mu', lower=0, upper=100, shape=_mu.shape)
    sigma = pm.Uniform('sigma', lower=0, upper=10, shape=_sigma.shape)

    # define likelihood
    y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=y)

    # inference
    start = pm.find_MAP()
    step = pm.Slice()
    trace = pm.sample(niter, step, start, random_seed=123, progressbar=True)
```

    Applied interval-transform to mu and added transformed mu_interval_ to model.
    Applied interval-transform to sigma and added transformed sigma_interval_ to model.
    100%|██████████| 1000/1000 [00:01<00:00, 918.42it/s]



```python
fig = plt.figure(figsize=(10,6))
ax1 = fig.add_subplot(1,2,1);
ax2 = fig.add_subplot(1,2,2);

ax1.hist(trace['mu'][-niter/2:,0], 25, histtype='step');
ax2.hist(trace['sigma'][-niter/2:,0], 25, histtype='step');
```

## Switchpoint analysis

This example comes from [Cam Davidson-Pilon's book](https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers)

We will be looking at count data---specifically at the frequency of text messages recieved over a period of time.  We will use the Poisson distribution to help use investigate this series of events.


```python
from IPython.display import Image
Image(filename='../poisson.png')
```


```python
fig = plt.figure(figsize=(12.5, 3.5))
ax = fig.add_subplot(111)
count_data = np.loadtxt("data/txtdata.csv")
n_count_data = len(count_data)
ax.bar(np.arange(n_count_data), count_data, color="#348ABD")
ax.set_xlabel("Time (days)")
ax.set_ylabel("count of text-msgs received")
ax.set_title("Did the user's texting habits change over time?")
ax.set_xlim(0, n_count_data);
```

These are time-course data.  Can we infer when a behavioral change occurred?

Lets denote the day with $i$ and text-message count by $C_i$,

$$ C_i \sim \text{Poisson}(\lambda) $$

   * Does the rate $\lambda$ change?
   * Is there a day (call it $\tau$) where the parameter $\lambda$ suddenly jumps to a higher value?
   * We are looking for a _switchpoint_ s.t.

$$
\lambda = 
\begin{cases}
\lambda_1  & \text{if } t \lt \tau \cr
\lambda_2 & \text{if } t \ge \tau
\end{cases}
$$

If no sudden change occurred and indeed $\lambda_1 = \lambda_2$, then the $\lambda$s posterior distributions should look about equal.

We are interested in inferring the unknown $\lambda$s. We are using Bayesian inference so we need priors. The *exponential* distribution provides a continuous density function for positive numbers.

\begin{align}
&\lambda_1 \sim \text{Exp}( \alpha ) \\\
&\lambda_2 \sim \text{Exp}( \alpha )
\end{align}

$\alpha$ is a *hyper-parameter*. In literal terms, it is a parameter that influences other parameters.

$$\frac{1}{N}\sum_{i=0}^N \;C_i \approx E[\; \lambda \; |\; \alpha ] = \frac{1}{\alpha}$$ 

What about $\tau$? Because of the noisiness of the data, it's difficult to pick out a priori when $\tau$ might have occurred. Instead, we can assign a *uniform prior belief* to every possible day. This is equivalent to saying

\begin{align}
& \tau \sim \text{DiscreteUniform(1,70) }\\\\
& \Rightarrow P( \tau = k ) = \frac{1}{70}
\end{align}

So what does this look like in PyMC3?


```python
import pymc3 as pm
import theano.tensor as tt

## assign lambdas and tau to stochastic variables
with pm.Model() as model:
    alpha = 1.0/count_data.mean()  # Recall count_data is the
                                   # variable that holds our txt counts
    lambda_1 = pm.Exponential("lambda_1", alpha)
    lambda_2 = pm.Exponential("lambda_2", alpha)
    tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)

## create a combined function for lambda (it is still a RV)    
with model:
    idx = np.arange(n_count_data) # Index
    lambda_ = pm.math.switch(tau >= idx, lambda_1, lambda_2)

## combine the data with our proposed data generation scheme    
with model:
    observation = pm.Poisson("obs", lambda_, observed=count_data)

## inference
with model:
    step = pm.Metropolis()
    trace = pm.sample(10000, tune=5000,step=step)
```

    Applied log-transform to lambda_1 and added transformed lambda_1_log_ to model.
    Applied log-transform to lambda_2 and added transformed lambda_2_log_ to model.
    100%|██████████| 10000/10000 [00:03<00:00, 3171.19it/s]



```python
## get the variables we want to plot from our trace
lambda_1_samples = trace['lambda_1']
lambda_2_samples = trace['lambda_2']
tau_samples = trace['tau']

# draw histogram of the samples:
fig = plt.figure(figsize=(12.5,10))
ax1 = fig.add_subplot(311)
ax2 = fig.add_subplot(312)
ax3 = fig.add_subplot(313)

for ax in [ax1,ax2]:
    ax.set_autoscaley_on(False)

## axis 1    
ax1.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85,
         label="posterior of $\lambda_1$", color="#A60628", normed=True)
ax1.legend(loc="upper left")
ax1.set_title(r"""Posterior distributions of the variables
    $\lambda_1,\;\lambda_2,\;\tau$""")
ax1.set_xlim([15, 30])
ax1.set_xlabel("$\lambda_1$ value")

## axis 2
ax2.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85,
         label="posterior of $\lambda_2$", color="#7A68A6", normed=True)
ax2.legend(loc="upper left")
ax2.set_xlim([15, 30])
ax2.set_xlabel("$\lambda_2$ value")

## axis 3
w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)
ax3.hist(tau_samples, bins=n_count_data, alpha=1,
         label=r"posterior of $\tau$",
         color="#467821", weights=w, rwidth=2.)
ax3.set_xticks(np.arange(n_count_data))

ax3.legend(loc="upper left")
ax3.set_ylim([0, .75])
ax3.set_xlim([35, len(count_data)-20])
ax3.set_xlabel(r"$\tau$ (in days)")
ax3.set_ylabel("probability");
plt.subplots_adjust(hspace=0.4)
```

### How do we plot the expection distribution of this model?


```python
fig = plt.figure(figsize=(12.5,5))
ax = fig.add_subplot(111)

N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_count_data)
for day in range(0, n_count_data):
    ix = day < tau_samples
    expected_texts_per_day[day] = (lambda_1_samples[ix].sum()
                                   + lambda_2_samples[~ix].sum()) / N

ax.plot(range(n_count_data), expected_texts_per_day, lw=4, color="#E24A33",
         label="expected number of text-messages received")
ax.set_xlim(0, n_count_data)
ax.set_xlabel("Day")
ax.set_ylabel("Expected # text-messages")
ax.set_title("Expected number of text-messages received")
ax.set_ylim(0, 60)
ax.bar(np.arange(len(count_data)), count_data, color="#348ABD", alpha=0.65,label="observed texts per day")

ax.legend(loc="upper left");
```

A distribution allows us to see the uncertainty in our estimates.  The two distribution for $\lambda$ are distinct.  Near day 45, there was a 50% chance that the user's behavior changed.

> In fact, the 45th day corresponds to Christmas, and I moved away to Toronto the next month, leaving a girlfriend behind

      -Cam Davidson Pilon



```python

```

9facdfaecb72b06f9eb337c0d05c25462b374bfd
ipynb
Jupyter Notebook
notebooks/introduction.ipynb
GalvanizeDataScience/probabilistic-programming-intro
9e007b6ae8e9d0ca2aff6c3c5e00c0aa4c012595
BSD-3-Clause
2017-02-08T04:18:46.000Z
2018-09-10T14:33:13.000Z
notebooks/introduction.ipynb
zipfian/probabilistic-programming-intro
9e007b6ae8e9d0ca2aff6c3c5e00c0aa4c012595
BSD-3-Clause
notebooks/introduction.ipynb
zipfian/probabilistic-programming-intro
9e007b6ae8e9d0ca2aff6c3c5e00c0aa4c012595
BSD-3-Clause
2017-02-09T01:06:04.000Z
2018-02-09T06:09:29.000Z
Qwen/Qwen-72B
1. YES
2. YES