<jupyter_start><jupyter_text>Bayesian Examples in Python<jupyter_code>import pandas as pd
import numpy as np

# Visualizations
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('fivethirtyeight')
plt.rcParams['font.size'] = 22
%matplotlib inline



import warnings
warnings.filterwarnings('ignore', category=FutureWarning)

import pymc3 as pm

# Helper functions
from utils import draw_pdf_contours, Dirichlet, plot_points, annotate_plot, add_legend, display_probs
# observations
animals = ['lions', 'tigers', 'bears']
c = np.array([5, 3, 1])



# hyperparameters (initially all equal)
alphas = np.array([1, 1, 1])<jupyter_output><empty_output><jupyter_text>[Calculate Expected Value](https://en.wikipedia.org/wiki/Categorical_distributionBayesian_inference_using_conjugate_prior)<jupyter_code>display_probs(dict(zip(animals, (alphas + c) / (c.sum() + alphas.sum()))))<jupyter_output>Species: lions    Prevalence: 50.00%.
Species: tigers   Prevalence: 33.33%.
Species: bears    Prevalence: 16.67%.<jupyter_text>Maximum A Posterior EstimationThe maximum a posterior (MAP) estimate is simply going to be the prevalence seen in the data. This is a frequentist way of viewing the world. $${\displaystyle \operatorname {arg\,max} \limits _{\mathbf {p} }p(\mathbf {p} \mid \mathbb {X} )={\frac {\alpha _{i}+c_{i}-1}{\sum _{i}(\alpha _{i}+c_{i}-1)}},\qquad \forall i\;\alpha _{i}+c_{i}>1}$$<jupyter_code>display_probs(dict(zip(animals, (alphas + c - 1) / sum(alphas + c - 1))))<jupyter_output>Species: lions    Prevalence: 55.56%.
Species: tigers   Prevalence: 33.33%.
Species: bears    Prevalence: 11.11%.<jupyter_text>Bayesian Model<jupyter_code>## Change n to match total prior observations

with pm.Model() as model:
    # Dirichlet hyperparameters are uniform hyperpriors
    hyperpriors = pm.Uniform('hyperpriors', shape = 3, observed = alphas)
    # Probabilities for each species
    parameters = pm.Dirichlet('parameters', a=hyperpriors, shape=3)
    # Observed data is a multinomial distribution with 6 trials
    observed_data = pm.Multinomial(
        'observed_data', n=9, p=parameters, shape=3, observed=c)
model
with model:
    # Sample from the posterior
    trace = pm.sample(draws=1000, chains=2, tune=500, discard_tuned_samples=False)
summary = pm.summary(trace)
summary.index = animals
summary
# Tuning samples
tune_df = pd.DataFrame(trace['parameters'][:1000], columns = animals)
tune_df['tune'] = True

# Samples after tuning
trace_df = pd.DataFrame(trace['parameters'][1000:], columns = animals)
trace_df['tune'] = False

all_df = pd.concat([tune_df, trace_df])

pvals = trace_df.iloc[:, :3].mean(axis = 0)
trace_df.head()
trace_df.iloc[:, :3].apply(lambda x: np.percentile(x, 5)).to_frame().rename(columns = {0: '5th percentile'})
## 95% range of uncertainity (Large amount of uncertinity given the small sampling size)
trace_df.iloc[:, :3].apply(lambda x: np.percentile(x, 95)).to_frame().rename(columns = {0: '95th percentile'})<jupyter_output><empty_output><jupyter_text>Diagonistic Plots<jupyter_code>ax = pm.plot_posterior(trace, varnames = ['parameters'], 
                       figsize = (20, 10), edgecolor = 'k');

plt.rcParams['font.size'] = 22
for i, a in enumerate(animals):
    ax[i].set_title(a);<jupyter_output><empty_output><jupyter_text>The traceplot shows a kernel density estimate and all the samples that were drawn on the right. We collapse the chains on th plots (combined = True) but in reality we drew 2 independent chains.<jupyter_code>prop_cycle = plt.rcParams['axes.prop_cycle']
cs = [x['color'] for x in list(prop_cycle)]

ax = pm.traceplot(trace, varnames = ['parameters'], figsize = (20, 10), combined = True,skip_first = 1000);
ax[0][0].set_title('KDE of Posteriors'); ax[0][1].set_title('Values in Trace');
add_legend(ax[0][0])
add_legend(ax[0][1])<jupyter_output><empty_output><jupyter_text>Maximum A Posteriori Result with PyMC3The MAP estimates are exactly the same as the observations. These are also the results that a frequentist would come up with!<jupyter_code>with model:
    # Find the maximum a posteriori estimate
    map_ = pm.find_MAP()
    
display_probs(dict(zip(animals, map_['parameters'])))
with model:
    samples = pm.sample_ppc(trace, samples = 1000)
    
dict(zip(animals, samples['observed_data'].mean(axis = 0)))
sample_df = pd.DataFrame(samples['observed_data'], columns = animals)

plt.figure(figsize = (22, 8))
for i, animal in enumerate(sample_df):
    plt.subplot(1, 3, i+1)
    sample_df[animal].value_counts().sort_index().plot.bar(color = 'r');
    plt.xticks(range(7), range(7), rotation = 0);
    plt.xlabel('Number of Times Seen'); plt.ylabel('Occurences');
    plt.title(f'1000 Samples for {animal}');<jupyter_output><empty_output><jupyter_text>Dirichlet Distributions<jupyter_code>draw_pdf_contours(Dirichlet(alphas))
annotate_plot()
draw_pdf_contours(Dirichlet(9*pvals))
annotate_plot()<jupyter_output><empty_output><jupyter_text>Probability of Next Observation<jupyter_code>next_obs = np.random.multinomial(n = 1, pvals = pvals, size = 10000)

# Data manipulation
next_obs = pd.melt(pd.DataFrame(next_obs, columns = ['Lions', 'Tigers', 'Bears'])).\
            groupby('variable')['value'].\
            value_counts(normalize=True).to_frame().\
             rename(columns = {'value': 'total'}).reset_index()
next_obs = next_obs.loc[next_obs['value'] == 1]

# Bar plot
next_obs.set_index('variable')['total'].plot.bar(figsize = (8, 6));
plt.title('Next Observation Likelihood');
plt.ylabel('Likelihood'); plt.xlabel('');
next_obs
all_df = pd.melt(all_df, id_vars=['tune'], 
                   var_name='animal', value_name='posterior')

plt.rcParams['font.size'] = 20
g = sns.FacetGrid(data = all_df, hue = 'tune', col = 'animal', size = 6)
g.map(sns.kdeplot, 'posterior');
l = plt.legend(prop={'size': 20});
l.set_title('Tune');
g.set_ylabels(label = 'Density')
summary<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup