---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.10.3
kernelspec:
  display_name: Python 3
  language: python
  name: pymc
---

```{code-cell} ipython3
import numpy as np
from scipy.stats import norm,uniform,lognorm
import matplotlib.pyplot as plt
import seaborn as sbn
import pandas as pd

np.random.seed(348348756)

plt.style.use('ggplot')
```

# Inference and Monte Carlo Markov Chains

Frequentists are used to the idea of hypothesis tests like 

\begin{eqnarray}
H_0: \mu = 0 \nonumber\\
H_1: \mu \neq 0 \nonumber
\end{eqnarray}

To proceed, we examine in expecation the moments of a parameter of
interest.  For example the variance of a mean $\bar{x}$ is computed as

$$
Var(x) = E_x\left[\left(x - E(x)\right)^2\right]
$$

an interesting take on this expecation is that it is with respect to
"data", since the true but unknown mean of $x$ exists and can be
written as $\mu_x$, making the above

$$
Var(x) = E_x\left[\left(x - \mu_x \right)^2\right]
$$

This variance can then be estimated from a sample having $N$
observations of data using

$$
\hat{Var(x)} = \frac{\sum_{i=1}^N(x_i-\bar{x})^2}{N}
$$

where $\bar{x} = \frac{\sum_{i=1} x_i }{N}$.

From there, confidence intervals, z scores, and other measures can be
computed inference.

The frequentist approach typically follows these steps:

1. Collect data
2. Define a model
3. Use data and model to calculate parameters of interest
4. Use parameters to test hypotheses

The Bayesian approach differs a bit

1. Collect data
2. Define a model
3. Use data and model to sample many times from the posterior distribution of parameters of interest (can be multivariate [many parameters])
4. Use the parameter samples to test hypothesis

It is really in steps (3) and (4) that things differ.  This notebook shows how this is implemented.

```{note}
Once can say that inference for a Bayesian is "attaching probabilities to statements".  These statements can be problem driven and need to be as formally defined as the hypothesis test above.
```

## A simple model

For having something concrete to play with, let's resurrect our simple example yet again:

```{code-cell} ipython3
data = norm(10,3).rvs(10)
sigma = 3          # Note this is the std of x
mu_prior = 8
sigma_prior = 1.5  # Note this is our prior on the std of mu

plt.hist(data,bins=20, density=True)
plt.show()
```

and we will again use our simple sampler:

```{code-cell} ipython3
# a fairly basic mh mcmc random walk sampler:
def sampler(data, samples=4, mu_init=9, sigma= 1, proposal_width=3, 
            mu_prior_mu=10, mu_prior_sd=1):
    mu_current = mu_init
    posterior = [mu_current]
    for i in range(samples):
        # suggest new position
        mu_proposal = norm(mu_current, proposal_width).rvs()

        # Compute likelihood by multiplying probabilities of each data point
        likelihood_current = norm(mu_current, sigma).pdf(data).prod()
        likelihood_proposal = norm(mu_proposal, sigma).pdf(data).prod()
        
        # Compute prior probability of current and proposed mu        
        prior_current = norm(mu_prior_mu, mu_prior_sd).pdf(mu_current)
        prior_proposal = norm(mu_prior_mu, mu_prior_sd).pdf(mu_proposal)
        
        p_current = likelihood_current * prior_current
        p_proposal = likelihood_proposal * prior_proposal
        
        # Accept proposal?
        c = np.min((1,p_proposal / p_current))
        
        accept = np.random.rand() <= c
                
        if accept:
            # Update position
            mu_current = mu_proposal
        
        posterior.append(mu_current)
        
    return np.array(posterior)
```

With that defined, let's take 5000 samples from the posterior:

```{code-cell} ipython3
chain = sampler(data,samples=5000,mu_init=9,sigma=sigma,proposal_width=3,
                mu_prior_mu=mu_prior,mu_prior_sd = sigma_prior)
```
 
### Using our chain for inference

In what follows, I assume we have already assessed chain convergence
and are satisfied with a burn-in of 100.

In the code below, I ask a number of questions (you might call this
inference) using our sampled posterior values, starting with ones most
similar to frequentist statistics:

```{code-cell} ipython3
chain_s = chain[100:]
```

Let's plot the histogram of our sample mean values:
```{code-cell} ipython3
plt.figure(figsize=(10,8))
plt.hist(chain_s, bins=75, density=True)
plt.show()
```

Here are some common inference-related questions a frequentist might ask:
```{code-cell}
print( "The mean is: ", chain_s.mean())
print( "The standard deviation is: ", chain_s.std())
print( "The median is: ", np.percentile(chain_s,50))
print( "The 95% Credible Interval is: ",np.percentile(chain_s,[2.5,97.5]))
print( "The probability that the mean is greater than 10: ", \
            chain_s[chain_s>10].shape[0]/chain_s.shape[0])
```

Some observations:

1. Notice that the MCMC analysis above provides thousands of potential
   mean estimates.  Some of these are closer to the fat part of the
   posterior while others are not.  We take an average of each of
   these potential mean values to estimate the mean of the posterior.
   We can also use these potential mean values to compute other
   moments of interest (like standard deviation).
2. We can also use the sample of means to compute non-parametric
   credibility intervals.  Note this doesn't make any distributional
   assumptions over the sample parameter values. For a Bayesian, the
   credible interval can be interpreted as the probability that the
   parameter of interest resides in the range given by the interval.
   For a frequentist confidence interval, it must be interpreted
   differently:

   ```{admonition} Frequentist Confidence Interval 
   Suppose we have a 95% confidence interval defined over [2.1, 5.9].  It is not correct to say that there is a 95% chance that the "true" parameter resides within 2.1 and 5.9.  Rather one needs to phrase it as follows:  If you collected many different data sets and computed your statistic, 95% of these computed statistics would reside in the range 2.1 and 5.9.  It is more of a statement about your confidence in the method you are using.  Recall, that for a frequentist, the unknown parameter is non-random so that either the parameter resides within 2.1 and 5.9 or not.
   ```
3. For "one-sided" hypothesis testing, we can simply compute the
   fraction of cases above the critical value that is being tested.

In short, samples from the posterior can provide nearly **any** type
of inference related question we might want to ask.  Overlaying the
credible interval on our histogram of sampled mean values, shows that
sample itself reveals what is and isn't most likely:

```{code-cell} ipython3
plt.figure(figsize=(10,8))
plt.hist(chain_s, density=True, bins=75)
plt.axvline(np.percentile(chain_s,[2.5,97.5])[0], c='k', ls='--')
plt.axvline(np.percentile(chain_s,[2.5,97.5])[1], c='k', ls='--');
```

While the last question isn't one we typically ask in a frequentist setting, we could by imposing normality from our model estimates and calculating areas under the pdf.  Alternatively, we could perform a bootstrapping excercise as we did earlier in the class.  The posterior pdf pictured above does not assume normality of our parameters and in that way is similar to the bootstrap approach for inference. 

**The Credible Interval** above can be interpreted as
> There is a 95% chance that $\theta$ is in the interval [8.84,9.95]

Note: the numbers in brackets will change depending on the credible interval printed above.

#### Classical Hypothesis testing

In frequentist statistics, we often perform test as follows:

$$
H_0 : \theta = 0 \\
H_1 : \theta \ne 0
$$

This isn't very easy to do in a Bayesian Context.  We can see that **no values** of our chain is exactly equal to 0 and this would likely be the case even if our mean for this problem was centered on zero. 

```{code-cell} ipython3
print( "Number of cases where theta=10: ", chain_s[chain_s==10].shape[0])
```

So does this mean we have 100% evidence that theta isn't 10? A
Bayesian would typically recast this question into an interval:

$$
H_0 : \theta \in 10 \pm .1 \\
H_1 : \theta \notin 10 \pm .1
$$

```{code-cell} ipython3
print( "Probability of H_0: ",chain_s[(chain_s<=10.1) & (chain_s>=9.9)].shape[0] / chain_s.shape[0])
```

So Bayesians like to attach probabilities to statements.  Additionally
note that due to our very small sample size, the prior is definitely
affecting our inference here.

#### Highest Posterior Density (HPD) Intervals

For skewed or multimodal posteriors, credible intervals may not be as
enlightening as they could be.  Suppose we have a MH markov chain and
our posterior pdf is plotted below.  The HPD interval, sets a minimum
probability threshold for parameters (e.g. posterior value must be
greater than .025) and examines the interval of values of $\theta$
that have more than a .025% chance of being observed.

```{note}
Skewed of multimodel samples from a posterior can be a sign that your chain hasn't converged or some other underlying pathology with your model.
```

Consider this skewed chain and an exercise of sliding a horizontal dotted black line down until it's y-value is
equal to the cutoff you as an analyst care about (e.g. .025):

```{code-cell} ipython3
skewed_chain = lognorm(1,8).rvs(10000)

plt.figure(figsize=(10,8))
plt.hist(skewed_chain,bins=100, density=True)
plt.axhline(.025, c='k', linestyle='--')
plt.axvline(12.74, c='k', linestyle='--')
plt.xlabel('$\\theta$')
plt.xlim((7,20))
plt.title('Posterior')
plt.show()
```

Note on the left hand size of the posterior, there is no region where the occurence of $\theta$ is less than .025.

```{code-cell} ipython3
print("The mean is: ", skewed_chain.mean())
print("The standard deviation is: ", skewed_chain.std())
print("The median is: ", np.percentile(skewed_chain,50))
print("The 95% Credible Interval is: ",np.percentile(skewed_chain,[2.5,97.5]))
```

Note the credible interval here is a bit strange, as there is large probability mass, near the minimum observed value. So turning to the HPD interval:

```{code-cell} ipython3
print("The highest posterior density interval is (approximately): ")
print("["+str(skewed_chain.min())+", 12.74")
print("Note: Min Chain Value is: ",skewed_chain.min())
```

Consider another case, a bimodal chain:

```{code-cell} ipython3
# suppose our chain was this
bimodal_chain = np.append(norm(2,.8).rvs(5000),norm(8,.5).rvs(5000))
plt.figure(figsize=(12,8))
plt.hist(bimodal_chain,bins=100, density=True)
plt.xlabel('$\\theta$')
plt.title('Posterior')
plt.show()
```

Applying the same principle, we can slide the horizontal dotted line
down until we find the transition from low ($<.025$) and high regions
of sampled parameters:

```{code-cell} ipython3
plt.figure(figsize=(12,7))
bimodal_chain = np.append(norm(2,.8).rvs(5000),norm(8,.5).rvs(5000))
plt.hist(bimodal_chain,bins=100, density=True)
plt.xlabel('$\\theta$')
plt.title('Posterior')
plt.axhline(.025/2,c='k',linestyle=':',label="95% HPD Cutoff for Bi-Modal Chain")
plt.axvline(bimodal_chain.mean(),c='r',label="Mean")
plt.axvline(np.percentile(bimodal_chain,2.5),c='g',linestyle='--',alpha=.5,label="95% CI")
plt.axvline(np.percentile(bimodal_chain,97.5),c='g',linestyle='--',alpha=.5)

plt.axvline(.05,c='k',linestyle='--',label="Approximate 95% HPD CI")
plt.axvline(4.05,c='k',linestyle='--')
plt.axvline(6.7,c='k',linestyle='--')
plt.axvline(9.39,c='k',linestyle='--')

plt.legend(loc='upper left', fontsize=14)
plt.show()
```

### The Chain and MAP

How do the moments of our chain (e.g. mean and median) relate to MAP (parameter value that maximes the posterior)?  To answer that question, we need to solve for MAP given the data:

```{code-cell} ipython3
def log_posterior(mu,data,sigma,mu_mu_prior,mu_sigma_prior):
    ln_prior = norm(mu_mu_prior,mu_sigma_prior).logpdf(mu)
    ln_like = norm(mu,sigma).logpdf(data).sum()
    return -1*(ln_like + ln_prior)
```

Using the negative log posterior, minimize to find MAP:
```{code-cell} ipython3
startvals = 9
from scipy.optimize import minimize
res = minimize(log_posterior,startvals, 
               args=(data,sigma,mu_prior,sigma_prior),options={'disp': True},
               method='BFGS')
```

And here are the results:
```{code-cell} ipython3
print("The MAP is: ", res.x[0])
print("The mean of our sampled values is: ", chain[100:].mean())
print("The median is: ", np.percentile(chain[100:],50))
```

For well behaved symmetric posteriors, the mean of the posterior will
usually be very close to MAP, but this doesn't necessarily have to be
the case.  Note that in the preceding figure (the bimodal posterior),
MAP is approximately 8 whereas the mean is approximately 5.