### Using PyMC3 for Bayesian Analysis
I will use PyMC3 to do some Bayesian analysis. I will follow the book by Davidson-Pilon.


```python
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
%matplotlib inline
```


```python
txt_data = np.loadtxt('data/txtdata.csv')
```


```python
n_txt_data = len(txt_data)
days = np.arange(n_txt_data)
```


```python
fig, ax = plt.subplots(1, 1, figsize = (25,9))
ax.bar(days, txt_data, color = 'xkcd:violet')
ax.set_xlabel(r'No. of Days', fontsize = 20)
ax.set_ylabel('No. of text messages', fontsize = 20)
ax.tick_params(axis = 'both', labelsize = 20 )
ax.set_facecolor('0.80')
ax.grid(which = 'both', axis = 'both', color = 'xkcd:white')
ax.set_axisbelow(True)

```

We want to model the texting behaviour of the author over various days. Since we are counting text messages and appropriate distribution function to consider would be the Poisson distribution. However since we want to check for texting behaviour we might want to consider two different Poisson distributions to get a better understanding.

Let us denote the text message count on the $i^{th}$ day as $C_i$. Since we are couunting here we can furthermore write:

$C_i = Poisson\left(\lambda\right)$

If we assume that $\lambda$ might change after somedayy, then we can assume two values $\lambda_1$ and say after $\tau$ days we have $\lambda_2$. The specific time $\tau$ is known as the *switchpoint*. 

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

Now in order to estimate $\lambda_1$ and $\lambda_2$ using Bayesian inference we need to have a prior distribution. We know that we can use an exponential distribution function to describe non-discrete positive numbers. This exponential distribution function will have a parameter of its own. let us call this parameter $\alpha$

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


Since the data is noisy it would be hard to decide a prior for $\tau$. We can thus consider a uniform prior in this case:

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


#### Using PyMC3


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

    WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.
    C:\Users\nathd\Anaconda3\lib\site-packages\h5py\__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
      from ._conv import register_converters as _register_converters



```python
with pm.Model() as model:
    alpha = 1.0/txt_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_txt_data - 1)
```

We can now create a function lambda_ which will allow us to switch between $\lambda_1$ and $\lambda_2$ after time $\tau$.


```python
with model:
    idx = np.arange(n_txt_data) # Index
    lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)
```

All the parameters $\lambda_1, \lambda_2, \text{ and } \tau$ are random number and so is lambda_ .


```python
with model:
    observation = pm.Poisson("obs", lambda_, observed=txt_data)
```

The variable `observation` combines our data, `count_data`, with our proposed data-generation scheme, given by the variable `lambda_`, through the `observed` keyword. 

The code below will be explained in Chapter 3, but I show it here so you can see where our results come from. One can think of it as a *learning* step. The machinery being employed is called *Markov Chain Monte Carlo* (MCMC), which I also delay explaining until Chapter 3. This technique returns thousands of random variables from the posterior distributions of $\lambda_1, \lambda_2$ and $\tau$. We can plot a histogram of the random variables to see what the posterior distributions look like. Below, we collect the samples (called *traces* in the MCMC literature) into histograms.


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

    Multiprocess sampling (2 chains in 2 jobs)
    CompoundStep
    >Metropolis: [tau]
    >Metropolis: [lambda_2_log__]
    >Metropolis: [lambda_1_log__]
    The number of effective samples is smaller than 25% for some parameters.



```python
lambda_1_samples = trace['lambda_1']
lambda_2_samples = trace['lambda_2']
tau_samples = trace['tau']
```


```python
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize = (15, 15))
ax1.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.6,
         label="posterior of $\lambda_1$", color="#A60628", normed=True)
ax1.legend(loc = 'best')
ax1.set_xlim([15, 30])
ax1.set_xlabel("$\lambda_1$")
ax1.set_facecolor('0.80')
ax1.grid(which = 'both', axis = 'both', color = 'xkcd:white', linestyle='dashed')
ax1.set_axisbelow(True)
ax2.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.5,
         label="posterior of $\lambda_2$", color="#7A68A6", normed=True)
ax2.legend(loc='best')
ax2.set_xlim([15, 30])
ax2.set_xlabel("$\lambda_2$")
ax2.set_facecolor('0.80')
ax2.grid(which = 'both', axis = 'both', color = 'xkcd:white', linestyle='dashed')
ax2.set_axisbelow(True)
w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)
ax3.hist(tau_samples, bins=n_txt_data, alpha=1,
         label=r"posterior of $\tau$",
         color="#467821", weights=w, rwidth=2.)
ax3.legend(loc='best')
ax3.set_xlim([35, len(txt_data)-20])
ax3.set_xlabel("Days")
ax3.set_ylabel('Probability')
ax3.set_facecolor('0.80')
ax3.grid(which = 'both', axis = 'both', color = 'xkcd:white', linestyle='dashed')
ax3.set_axisbelow(True)
```

We see that the values of $\lambda_1$ and $\lambda_2$ is around 18 and 23 respectively. Additionally we see that the change in text messages received takes place over a span of 4 days, with the likelihood that the change in the number of text takes place on the 45th day is around 50 %. We would like to see how this change of number of text messages received takes place over time.


```python
fig, ax = plt.subplots(1,1, figsize = (15, 5))

# tau_samples, lambda_1_samples, lambda_2_samples contain
# N samples from the corresponding posterior distribution
N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_txt_data) 
for day in range(n_txt_data):
    #ix is the boolean index of all tay samples corresponding to the switch point occuring prior to the value of day
    ix = day<tau_samples
    # Each posterior sample corresponds to a value for tau.
    # for each day, that value of tau indicates whether we're "before"
    # (in the lambda1 "regime") or
    #  "after" (in the lambda2 "regime") the switchpoint.
    # by taking the posterior sample of lambda1/2 accordingly, we can average
    # over all samples to get an expected value for lambda on that day.
    # As explained, the "message count" random variable is Poisson distributed,
    # and therefore lambda (the poisson parameter) is the expected value of
    # "message count".
    expected_texts_per_day[day] = (lambda_1_samples[ix].sum()
                                   + lambda_2_samples[~ix].sum()) / N
ax.plot(range(n_txt_data), expected_texts_per_day, '-', color = 'xkcd:magenta', linewidth = 5, alpha = 0.7, label = 'Expected Text Per Day')
ax.bar(days, txt_data, color = 'xkcd:violet', label = 'Observed texts per day', alpha = 0.7)
ax.legend(loc = 'best')
ax.set_xlabel(r'Days', fontsize = 20)
ax.set_ylabel('No. of text messages', fontsize = 20)
ax.tick_params(axis = 'both', labelsize = 20 )
ax.set_facecolor('0.80')
ax.grid(which = 'both', axis = 'both', color = 'xkcd:white')
ax.set_axisbelow(True)
```


```python

```

e42246f8c2994f285656108f4fde3c79e7f9c4be
ipynb
Jupyter Notebook
Markov Chain Monte Carlo/Bayesian Statistics using PyMC3.ipynb
nathdipankar/Bayesian-Statistics
444d0476e61e2ce61ae39bf0bb30c6d0e1b3d2c9
MIT
Markov Chain Monte Carlo/Bayesian Statistics using PyMC3.ipynb
nathdipankar/Bayesian-Statistics
444d0476e61e2ce61ae39bf0bb30c6d0e1b3d2c9
MIT
Markov Chain Monte Carlo/Bayesian Statistics using PyMC3.ipynb
nathdipankar/Bayesian-Statistics
444d0476e61e2ce61ae39bf0bb30c6d0e1b3d2c9
MIT
Qwen/Qwen-72B
1. YES
2. YES