<jupyter_start><jupyter_text>Estimation and Comparison with Probabilistic Programming<jupyter_code>import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import pymc3 as pm
from ipywidgets import interact
import arviz as az
import pandas as pd
import janitor
from utils import ECDF
import holoviews as hv
import hvplot.pandas

hv.extension("bokeh")


%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

sns.set_style('white')
sns.set_context('talk')<jupyter_output><empty_output><jupyter_text>Learning Objectives of Part 3 onwards 1. Consolidate your knowledge of the Bayesian model building workflow and use probabilistic programming for parameter estimation;2. Use probabilistic programming for hypothesis testing.  Bayesian parameter estimation using PyMC3 Well done! You've learnt the basics of Bayesian model building. The steps are1. To completely specify the model in terms of _probability distributions_. This includes specifying     - what the form of the sampling distribution of the data is _and_     - what form describes our _uncertainty_ in the unknown parameters (This formulation is adapted from [Fonnesbeck's workshop](https://github.com/fonnesbeck/intro_stat_modeling_2017/blob/master/notebooks/2.%20Basic%20Bayesian%20Inference.ipynb) as Chris said it so well there).2. Calculate the _posterior distribution_.In the above, the form of the sampling distribution of the data was Binomial (described by the likelihood) and the uncertainty around the unknown parameter $p$ captured by the prior. Now it is time to do the same using the **probabilistic programming language** PyMC3. There's _loads of cool stuff_ about PyMC3 and this paradigm, two of which are- _probabililty distributions_ are first class citizens, in that we can assign them to variables and use them intuitively to mirror how we think about priors, likelihoods & posteriors.- PyMC3 calculates the posterior for us: this is fancy math done for lazy programmers!Under the hood, PyMC3 will compute the posterior using a sampling based approach called Markov Chain Monte Carlo (MCMC) or Variational Inference. Check the [PyMC3 docs](https://docs.pymc.io/) for more on these. From this notebook onwards, we have prepared a series of examples that will show you how to perform inference and prediction in a variety of problems. Hopefully, you will see that at the end of the day, everything we do boils down to estimation of some kind, and by doing it in a Bayesian setting, we avoid many pitfalls that come from blindly following canned statistical procedures.  Example: Click-Through Rates A common experiment in tech data science is to test a product change and see how it affects a metric that you're interested in. Say that we don't think enough people are clicking a button on my website, but we hypothesize that it's because the button is a similar color to the background of the page, meaning they're a difficult time finding the buttong to click.We can serve up two pages, and randomly send some customers to each: the first the original page, the second a page that is identical, except that it has a button that is of higher contrast and see if more people click through. This is commonly referred to as an A/B test and the metric of interest is click-through rate (CTR), what proportion of people click through. We will use this example to help us build familiarity with PyMC3 mechanics.  Load DataLet's first load some data.<jupyter_code>ctr = (
    pd.read_csv('../data/clickthrough.csv', index_col=0)
    .label_encode('group')  # FYI: this is a pyjanitor function    
)
ctr.sample(10)<jupyter_output><empty_output><jupyter_text>Build Model: Estimate $p$ for control groupNow it's time to build our probability model. Noticing that our model of having a constant click-through rate resulting in click or not is a biased coin flip:- the sampling distribution is binomial and we need to encode this in the likelihood;- there is a single parameter $p$ that we need to describe the uncertainty around.     - Its value must be bound between 0 and 1, so we must use a probability distribution that takes on these bounds.    - Having not seen the data, we must express what we believe about the parameter $p$ -- this is nothing more than assigning credibility points to the number line for the value $p$. As a starter, we will show you how to express **"equal credibility from 0 to 1"**, by using the Uniform distribution.  Model DefinitionThese are the ingredients for the model so let's now build it.<jupyter_code>control_df = ctr.query("group == 'control'")
with pm.Model() as model1_bernoulli:
    p = pm.Uniform("p", lower=0, upper=1)
    like = pm.Bernoulli("likelihood", p=p, observed=control_df['clicks'])<jupyter_output><empty_output><jupyter_text>There is an alternative way to build this model, taking advantage of the fact that the sum of Bernoulli-distributed data follows a Binomial distribution. The syntax would look like this below, annotated with other modifications to guide you along:<jupyter_code># Build model of p_a
with pm.Model() as model1_binomial:
    # Prior on p
    p = pm.Uniform('p')  # defaults for the uniform distribution are lower=0 and upper=1
    # Binomial Likelihood
    like = pm.Binomial(
        'likelihood', 
        n=len(control_df), 
        p=p, 
        observed=len(control_df.query("clicks == 1"))
    )<jupyter_output><empty_output><jupyter_text>Note how the data have to be in a different shape, though. With the Bernoulli likelihood, we need every single success/failure to be recorded. With the Binomial likelihood, we only need the summary statistics: number of trials, and number of successes.  Sample from PosteriorIt's now time to sample from the posterior using PyMC3. You'll also plot the posterior:<jupyter_code>with model1_bernoulli:
    samples_bernoulli = pm.sample(2000, tuning=1000)
with model1_binomial:
    samples_binomial = pm.sample(2000, tuning=1000)<jupyter_output><empty_output><jupyter_text>Model CheckingNow, let's use ArviZ to perform some visual diagnostics.<jupyter_code>import arviz as az

# Posterior plot for bernoulli model
az.plot_posterior(samples_bernoulli, kind='hist', round_to=2);
# Posterior plot for binomial model
az.plot_posterior(samples_binomial, kind='hist', round_to=2);<jupyter_output><empty_output><jupyter_text>Notice how we get the same results using the two model formulations (Bernoulli- vs. binomial-distributed likelihoods).  DiscussionInterpret the posterior ditribution. What would your tell the non-technical manager of your growth team about the CTR?  Build Model: Compare $p$ for control and test groupsHaving built the model for the control group, let's now extend it to compare the control and test groups.  Hands-on: Build ModelModify the first model (go ahead and copy/paste it here!) to estimate $p_{control}$ and $p_{test}$. Let's use Binomial likelihoods, for a bit of variety.**Hint:** You will probably want to have two different `p`s (e.g. `p_control` and `p_test`) as well as two likelihoods (`like_control`, and `like_test`)<jupyter_code># Create test_df, just as we created control_df above.
test_df = ctr.query("group == 'test'")

# Give your model a variable name
with pm.Model() as model2:
    # Copy p and likelihood for control group.
    p_control = pm.Uniform('p_control')
    like_control = pm.Binomial(
        'like_control', 
        n=len(control_df), 
        p=p_control, 
        observed=len(control_df.query("clicks == 1"))
    )
    
    # Modify the above p and likelihood for test group.
    p_test = pm.Uniform('p_test')
    like_test = pm.Binomial(
        'like_test', 
        n=len(test_df), 
        p=p_test, 
        observed=len(test_df.query("clicks == 1"))
    )

    # We will also explicitly compute the difference between
    # p_control and p_test. 
    # This shows you that we can do math on probability distributions!
    p_diff = pm.Deterministic("p_diff", p_test - p_control)<jupyter_output><empty_output><jupyter_text>Sample from Posterior<jupyter_code>with model2:
    trace = pm.sample(2000)<jupyter_output><empty_output><jupyter_text>Model Checking<jupyter_code>az.plot_posterior(trace, round_to=2);<jupyter_output><empty_output><jupyter_text>DiscussionHow did the test group compare to the control group?  Hypothesis Testing<NAME> wrote [a blog post](http://allendowney.blogspot.com/2016/06/there-is-still-only-one-test.html) that illustrates how all of statistical inference boils down into a single framework: test statistic, null model, and "probability of unexpectedness" (p-value) under that null model. All of the special statistical tests with fancy author names are merely particular cases of this framework.  Code-along: Translating probability distributions into things that matterWith Bayesian methods, we will go out on a limb to say that the p-value doesn't matter at all: posterior distributions are all we need, and we can use posterior distributions to do some cool stuff.Let's make use of the posterior distribution on differences to calculate how much money we expect to gain.If we know that our customers on average spend 25 USD after clicking, and 0 USD if they do not click, we can then simulate the distribution of expected increase of revenue over 1 million customers.<jupyter_code>x, y = ECDF(trace['p_diff'] * 25 * 1E6)
plt.plot(x, y)
plt.xlabel("Expected revenue increase")
plt.ylabel("Cumulative probability")
sns.despine()<jupyter_output><empty_output><jupyter_text>Example: Baseball PlayersWe're now going to switch to a different dataset: that of baseball players and their batting stats. The goal of this analysis is to identify which player we want to target to make an offer to join our baseball team. To simplify the problem, we are going to make a decision on the basis of just batting average and their salaries.  Learning Objectives1. Reinforce Binomial/bernoulli generative story.1. Illustrate how to use broadcasting to avoid writing for-loops. 1. How to write a hierarchical model.  Load DataLet's load baseball player data.<jupyter_code>%qtconsole
players = (
    pd.read_csv("../data/baseballdb/core/Batting.csv")
    .clean_names()
    .query("yearid == 2016")
    .select_columns(['playerid', 'ab', 'h'])
    .groupby('playerid').sum()
)

salaries = (
    pd.read_csv("../data/baseballdb/core/Salaries.csv")
    .clean_names()
    .query("yearid == 2016")
    .select_columns(["playerid", "salary"])
    .groupby('playerid').mean()
)

data = (
    players.join(salaries)
    .dropna()
    .reset_index().label_encode("playerid").set_index("playerid")
)
data.head(10)<jupyter_output><empty_output><jupyter_text>Code-along: Build Model while avoiding for-LoopsPreviously, we were able to compare two groups by copy/pasting code. That's not an ideal way to construct our model here. We're going to show you how you can avoid for-loops by taking advantage of broadcasting.  Model DefinitionOnce again, we see the Bernoulli distribution data generation story at play here: we have a number of trials (`ab`), with a number of successes (`h`), from which we want to estimate a player-specific property: `p`, the probability of hitting a pitch.We are going to introduce a new distribution here: the `beta` distribution, which has the same desirable properties as the Uniform, but provides richer information.<jupyter_code>with pm.Model() as pitch_model:
    # we construct a vector of beta distributions, one for each player.
    p = pm.Beta("p", alpha=1, beta=1, shape=(len(data)))
    like = pm.Binomial("like", p=p, n=data['ab'], observed=data['h'])
    
    # Let's also construct a "p per salary" metric, based on their salary.
    # We want highest p for lowest salary.
    pps = pm.Deterministic("pps", p / data['salary'])<jupyter_output><empty_output><jupyter_text>Sample from Posterior<jupyter_code># Use the correct model variable here
with pitch_model:
    # Give your trace a name
    trace_batting = pm.sample(2000)<jupyter_output><empty_output><jupyter_text>Model InterpretationBecause we have 805 players, we're going to create a custom visualization that lets us scrub through our players' posterior distributions.<jupyter_code>from ipywidgets import SelectMultiple, HBox, VBox, Select, fixed

def dict2tuple(d):
    return list(*zip(*zip(d.items())))

def inversedict(d):
    return {v:k for k, v in d.items()}

playerid_mapping = data['playerid_enc'].to_dict()
playerid_inverse = inversedict(playerid_mapping)

player_select = SelectMultiple(
    options=dict2tuple(playerid_mapping),
    value=(0,)
)

metric_select = Select(
    options=['p', 'pps'],
)
def scrub_posterior(player_encs: list, metric: str, trace):
    fig, ax = plt.subplots()
    for enc in player_encs:
        x, y = ECDF(trace[metric][:, enc])
        ax.plot(x, y, label=playerid_inverse[enc])
    ax.legend()
    ax.set_xlabel("p per salary")
    ax.set_ylabel('cumulative probability')
        
    sns.despine()
interact(scrub_posterior, player_encs=player_select, metric=metric_select, trace=fixed(trace_batting));<jupyter_output><empty_output><jupyter_text>Hands-on: Finding our ideal playerUsing the visualization, find a bunch of players whom you think would be great to scout, based on their PPS.  Hierarchical model**Discussion:** Is it reasonable that some players have a `p` (probability of hitting) value that ranges from as small as 0 to as high as 1? Under what assumption would this be reasonable? Under what assumption would this be *unreasonable*?  Code-along: Model DefinitionTo build a hierarchical model, we essentially need to define a prior on the parameters of the Beta distribution. These new priors now become "hyperpriors". How do we pick hyperpriors then? Some qualitative rules of thumb:1. Pick a distribution with the correct **support**. For example, if a parameter can/should never take negative values, then pick a distribution with only positive support.1. The distribution should be **weakly informative**. For example, rather than a uniform distribution hyperprior, one might choose a wide distribution (relative to the relevant scale of the parameters) that spans 1 order of magnitude difference, but still has uneven allocation of credibility. Use the PyMC3 distribution library to your advantage! There are pictures that show you the shapes of the distributions, which can be helpful in narrowing your choices.<jupyter_code>with pm.Model() as hierarchical_baseball:
    a_prior = pm.Exponential("a_prior", lam=1/29.)
    b_prior = pm.Exponential("b_prior", lam=1/29.)
    
    p = pm.Beta("p", alpha=a_prior, beta=b_prior, shape=(len(data)))
    
    like = pm.Binomial("like", p=p, n=data['ab'], observed = data['h'])
    
    pps = pm.Deterministic("pps", p / data['salary'])<jupyter_output><empty_output><jupyter_text>Sample from Posterior<jupyter_code>with hierarchical_baseball:
    trace_hierarchical = pm.sample(2000)<jupyter_output><empty_output><jupyter_text>While sampling happens, do you have any questions? Please feel free to ask.  Hands-on: Model Interpretation<jupyter_code>interact(scrub_posterior, player_encs=player_select, metric=metric_select, trace=fixed(trace_hierarchical));<jupyter_output><empty_output><jupyter_text>**Discussion:** 1. What do you notice about the posterior distributions on `p`?1. What do you notice about the posterior distributions on `pps`?  Model Building Workflow**Things we did:**1. Get data into correct shape.1. Build naive model, followed by more complex model.1. Build a simple loss function into our model.**Things we didn't do:**1. Posterior predictive checks.**Bonus things:**1. Build custom posterior explorers. You can use your favourite framework of choice.  Hands-on: Parameter estimation II -- the mean of a population In this exercise, you'll calculate the  posterior mean beak depth of Galapagos finches in a given species. First you'll load the data and subset wrt species:<jupyter_code># Import and view head of data
df_12 = pd.read_csv('../data/finch_beaks_2012.csv')
df_fortis = df_12.loc[df_12['species'] == 'fortis']
df_scandens = df_12.loc[df_12['species'] == 'scandens']<jupyter_output><empty_output><jupyter_text>To specify the full probabilty model, you need- a likelihood function for the data &- priors for all unknowns.What is the likelihood here? Let's plot the measurements below and see that they look approximately Gaussian/normal so you'll use a normal likelihood $y_i\sim \mathcal{N}(\mu, \sigma^2)$. The unknowns here are the mean $\mu$ and standard deviation $\sigma$ and we'll use weakly informative priors on both- a normal prior for $\mu$ with mean $10$ and standard deviation $5$;- a uniform prior for $\sigma$ bounded between $0$ and $10$.We can discuss biological reasons for these priors also but you can also test that the posteriors are relatively robust to the choice of prior here due to the amount of data.<jupyter_code>sns.distplot(df_fortis['blength']);
df_fortis.head()
with pm.Model() as model:
    # Prior for mean & standard deviation
    μ_1 = pm.Normal('μ_1', mu=10, sd=5)
    σ_1 = pm.Lognormal('σ_1', 0, 10)
    # Gaussian Likelihood
    y_1 = pm.Normal('y_1', mu=μ_1, sd=σ_1, observed=df_fortis['blength'])
# bust it out & sample
with model:
    samples = pm.sample(2000, njobs=1)
az.plot_posterior(samples);<jupyter_output><empty_output><jupyter_text>2. Bayesian Hypothesis testing  Bayesian Hypothesis testing I: A/B tests on click through rates Assume we have a website and want to redesign the layout (*A*) and test whether the new layout (*B*) results in a higher click through rate. When people come to our website we randomly show them layout *A* or *B* and see how many people click through for each. First let's generate the data we need:<jupyter_code># click-through rates
p_a = 0.15
p_b = 0.20
N = 1000
n_successes_a = np.sum(np.random.uniform(size=N) <= p_a)
n_successes_b = np.sum(np.random.uniform(size=N) <= p_b)<jupyter_output><empty_output><jupyter_text>Once again, we need to specify our models for $p_a$ and $p_b$. Each will be the same as the CTR example above- Binomial likelihoods- uniform priors on $p_a$ and $_p$.We also want to calculate the posterior of the difference $p_a-p_b$ and we do so using `pm.Deterministic()`, which specifies a deterministic random variable, i.e., one that is completely determined by the values it references, in the case $p_a$ & $p_b$.We'll now build the model:<jupyter_code>with pm.Model() as Model:
    # Prior on p
    prob_a = pm.Uniform('p_a')
    prob_b = pm.Uniform('p_b')
    # Binomial Likelihood
    y_a = pm.Binomial('y_a', n=N, p=prob_a, observed=n_successes_a)
    y_b = pm.Binomial('y_b', n=N, p=prob_b, observed=n_successes_b)
    diff_clicks = pm.Deterministic('diff_clicks', prob_a-prob_b)<jupyter_output><empty_output><jupyter_text>Sample from the posterior and plot them:<jupyter_code>with Model:
    samples = pm.sample(2000, njobs=1)
az.plot_posterior(samples, kind='hist');<jupyter_output><empty_output><jupyter_text>Hands-on: Bayesian Hypothesis testing II -- beak lengths difference between species **Task**: Determine whether the mean beak length of the Galapogas finches differs between species. For the mean of each species, use the same model as in previous hand-on section:- Gaussian likelihood;- Normal prior for the means;- Uniform prior for the variances.Also calculate the difference between the means and, for bonus points, the _effect size_, which is the difference between the means divided by the pooled standard deviations = $\sqrt{(\sigma_1^2+\sigma_2^2)/2}$. Hugo will talk through the importance of the _effect size_.Don't forget to sample from the posteriors and plot them!<jupyter_code>with pm.Model() as model:
    # Priors for means and variances
    μ_1 = pm.Normal('μ_1', mu=10, sd=5)
    σ_1 = pm.Uniform('σ_1', 0, 10)
    μ_2 = pm.Normal('μ_2', mu=10, sd=5)
    σ_2 = pm.Uniform('σ_2', 0, 10)
    # Gaussian Likelihoods
    y_1 = pm.Normal('y_1', mu=μ_1, sd=σ_1, observed=df_fortis['blength'])
    y_2 = pm.Normal('y_2', mu=μ_2, sd=σ_2, observed=df_scandens['blength'])
    # Calculate the effect size and its uncertainty.
    diff_means = pm.Deterministic('diff_means', μ_1 - μ_2)
    pooled_sd = pm.Deterministic('pooled_sd', 
                                 np.sqrt(np.power(σ_1, 2) + 
                                         np.power(σ_2, 2)) / 2)
    effect_size = pm.Deterministic('effect_size', 
                                   diff_means / pooled_sd)
# bust it out & sample
with model:
    samples = pm.sample(2000, njobs=1)
az.plot_posterior(samples, var_names=['μ_1', 'μ_2', 'diff_means', 'effect_size'], kind='hist');
# Estimation and Comparison with Probabilistic Programming

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import pymc3 as pm
from ipywidgets import interact
import arviz as az
import pandas as pd
import janitor
from utils import ECDF
import holoviews as hv
import hvplot.pandas

hv.extension("bokeh")


%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

sns.set_style('white')
sns.set_context('talk')

## Learning Objectives of Part 3 onwards

1. Consolidate your knowledge of the Bayesian model building workflow and use probabilistic programming for parameter estimation;
2. Use probabilistic programming for hypothesis testing.

## Bayesian parameter estimation using PyMC3

Well done! You've learnt the basics of Bayesian model building. The steps are
1. To completely specify the model in terms of _probability distributions_. This includes specifying 
    - what the form of the sampling distribution of the data is _and_ 
    - what form describes our _uncertainty_ in the unknown parameters (This formulation is adapted from [Fonnesbeck's workshop](https://github.com/fonnesbeck/intro_stat_modeling_2017/blob/master/notebooks/2.%20Basic%20Bayesian%20Inference.ipynb) as Chris said it so well there).
2. Calculate the _posterior distribution_.

In the above, the form of the sampling distribution of the data was Binomial (described by the likelihood) and the uncertainty around the unknown parameter $p$ captured by the prior.

Now it is time to do the same using the **probabilistic programming language** PyMC3. There's _loads of cool stuff_ about PyMC3 and this paradigm, two of which are
- _probabililty distributions_ are first class citizens, in that we can assign them to variables and use them intuitively to mirror how we think about priors, likelihoods & posteriors.
- PyMC3 calculates the posterior for us: this is fancy math done for lazy programmers!

Under the hood, PyMC3 will compute the posterior using a sampling based approach called Markov Chain Monte Carlo (MCMC) or Variational Inference. Check the [PyMC3 docs](https://docs.pymc.io/) for more on these. 

From this notebook onwards, we have prepared a series of examples that will show you how to perform inference and prediction in a variety of problems. Hopefully, you will see that at the end of the day, everything we do boils down to estimation of some kind, and by doing it in a Bayesian setting, we avoid many pitfalls that come from blindly following canned statistical procedures.

## Example: Click-Through Rates

A common experiment in tech data science is to test a product change and see how it affects a metric that you're interested in. 

Say that we don't think enough people are clicking a button on my website, but we hypothesize that it's because the button is a similar color to the background of the page, meaning they're a difficult time finding the buttong to click.

We can serve up two pages, and randomly send some customers to each: the first the original page, the second a page that is identical, except that it has a button that is of higher contrast and see if more people click through. 

This is commonly referred to as an A/B test and the metric of interest is click-through rate (CTR), what proportion of people click through. We will use this example to help us build familiarity with PyMC3 mechanics.

### Load Data

Let's first load some data.

ctr = (
    pd.read_csv('../data/clickthrough.csv', index_col=0)
    .label_encode('group')  # FYI: this is a pyjanitor function    
)
ctr.sample(10)

### Build Model: Estimate $p$ for control group

Now it's time to build our probability model. Noticing that our model of having a constant click-through rate resulting in click or not is a biased coin flip:

- the sampling distribution is binomial and we need to encode this in the likelihood;
- there is a single parameter $p$ that we need to describe the uncertainty around. 
    - Its value must be bound between 0 and 1, so we must use a probability distribution that takes on these bounds.
    - Having not seen the data, we must express what we believe about the parameter $p$ -- this is nothing more than assigning credibility points to the number line for the value $p$. As a starter, we will show you how to express **"equal credibility from 0 to 1"**, by using the Uniform distribution.

#### Model Definition

These are the ingredients for the model so let's now build it.

control_df = ctr.query("group == 'control'")

with pm.Model() as model1_bernoulli:
    p = pm.Uniform("p", lower=0, upper=1)
    like = pm.Bernoulli("likelihood", p=p, observed=control_df['clicks'])

There is an alternative way to build this model, taking advantage of the fact that the sum of Bernoulli-distributed data follows a Binomial distribution. The syntax would look like this below, annotated with other modifications to guide you along:

# Build model of p_a
with pm.Model() as model1_binomial:
    # Prior on p
    p = pm.Uniform('p')  # defaults for the uniform distribution are lower=0 and upper=1
    # Binomial Likelihood
    like = pm.Binomial(
        'likelihood', 
        n=len(control_df), 
        p=p, 
        observed=len(control_df.query("clicks == 1"))
    )

Note how the data have to be in a different shape, though. With the Bernoulli likelihood, we need every single success/failure to be recorded. With the Binomial likelihood, we only need the summary statistics: number of trials, and number of successes.

#### Sample from Posterior


It's now time to sample from the posterior using PyMC3. You'll also plot the posterior:

with model1_bernoulli:
    samples_bernoulli = pm.sample(2000, tuning=1000)

with model1_binomial:
    samples_binomial = pm.sample(2000, tuning=1000)

#### Model Checking

Now, let's use ArviZ to perform some visual diagnostics.

import arviz as az

# Posterior plot for bernoulli model
az.plot_posterior(samples_bernoulli, kind='hist');

# Posterior plot for binomial model
az.plot_posterior(samples_binomial, kind='hist');

Notice how we get the same results using the two model formulations (Bernoulli- vs. binomial-distributed likelihoods).

#### Discussion

Interpret the posterior ditribution. What would your tell the non-technical manager of your growth team about the CTR?

### Build Model: Compare $p$ for control and test groups

Having built the model for the control group, let's now extend it to compare the control and test groups.

#### Hands-on: Build Model

Modify the first model (go ahead and copy/paste it here!) to estimate $p_{control}$ and $p_{test}$. Let's use Binomial likelihoods, for a bit of variety.

**Hint:** You will probably want to have two different `p`s (e.g. `p_control` and `p_test`) as well as two likelihoods (`like_control`, and `like_test`)

# Create test_df, just as we created control_df above.
test_df = ctr.query("group == 'test'")

# Give your model a variable name
with pm.Model() as model2:
    # Copy p and likelihood for control group.
    p_control = pm.Uniform('p_control')
    like_control = pm.Binomial(
        'like_control', 
        n=len(control_df), 
        p=p_control, 
        observed=len(control_df.query("clicks == 1"))
    )
    
    # Modify the above p and likelihood for test group.
    p_test = pm.Uniform('p_test')
    like_test = pm.Binomial(
        'like_test', 
        n=len(test_df), 
        p=p_test, 
        observed=len(test_df.query("clicks == 1"))
    )

    # We will also explicitly compute the difference between
    # p_control and p_test. 
    # This shows you that we can do math on probability distributions!
    p_diff = pm.Deterministic("p_diff", p_test - p_control)

#### Sample from Posterior

with model2:
    trace = pm.sample(2000)

#### Model Checking

az.plot_posterior(trace, round_to=2);

#### Discussion

How did the test group compare to the control group?

#### Hypothesis Testing

<NAME> wrote [a blog post](http://allendowney.blogspot.com/2016/06/there-is-still-only-one-test.html) that illustrates how all of statistical inference boils down into a single framework: test statistic, null model, and "probability of unexpectedness" (p-value) under that null model. All of the special statistical tests with fancy author names are merely particular cases of this framework.

#### Code-along: Translating probability distributions into things that matter

With Bayesian methods, we will go out on a limb to say that the p-value doesn't matter at all: posterior distributions are all we need, and we can use posterior distributions to do some cool stuff.

Let's make use of the posterior distribution on differences to calculate how much money we expect to gain.

If we know that our customers on average spend 25 USD after clicking, and 0 USD if they do not click, we can then simulate the distribution of expected increase of revenue over 1 million customers.

x, y = ECDF(trace['p_diff'] * 25 * 1E6)
plt.plot(x, y)
plt.xlabel("Expected revenue increase")
plt.ylabel("Cumulative probability")
sns.despine()

## Example: Baseball Players

We're now going to switch to a different dataset: that of baseball players and their batting stats. The goal of this analysis is to identify which player we want to target to make an offer to join our baseball team. 

To simplify the problem, we are going to make a decision on the basis of just batting average and their salaries.

### Learning Objectives

1. Reinforce Binomial/bernoulli generative story.
1. Illustrate how to use broadcasting to avoid writing for-loops. 
1. How to write a hierarchical model.

### Load Data

Let's load baseball player data.

players = (
    pd.read_csv("../data/baseballdb/core/Batting.csv")
    .clean_names()
    .query("yearid == 2016")
    .select_columns(['playerid', 'ab', 'h'])
    .groupby('playerid').sum()
)

salaries = (
    pd.read_csv("../data/baseballdb/core/Salaries.csv")
    .clean_names()
    .query("yearid == 2016")
    .select_columns(["playerid", "salary"])
    .groupby('playerid').mean()
)

data = (
    players.join(salaries)
    .dropna()
    .reset_index().label_encode("playerid").set_index("playerid")
)
data

### Code-along: Build Model while avoiding for-Loops

Previously, we were able to compare two groups by copy/pasting code. That's not an ideal way to construct our model here. We're going to show you how you can avoid for-loops by taking advantage of broadcasting.

#### Model Definition

Once again, we see the Bernoulli distribution data generation story at play here: we have a number of trials (`ab`), with a number of successes (`h`), from which we want to estimate a player-specific property: `p`, the probability of hitting a pitch.

We are going to introduce a new distribution here: the `beta` distribution, which has the same desirable properties as the Uniform, but provides richer information.

with pm.Model() as pitch_model:
    # we construct a vector of beta distributions, one for each player.
    p = pm.Beta("p", alpha=1, beta=1, shape=(len(data)))
    like = pm.Binomial("like", p=p, n=data['ab'], observed = data['h'])
    
    # Let's also construct a "p per salary" metric, based on their salary.
    # We want highest p for lowest salary.
    pps = pm.Deterministic("pps", p / data['salary'])
   

#### Sample from Posterior

# Use the correct model variable here
with pitch_model:
    # Give your trace a name
    trace_batting = pm.sample(2000)

#### Model Interpretation

Because we have 805 players, we're going to create a custom visualization that lets us scrub through our players' posterior distributions.

from ipywidgets import SelectMultiple, HBox, VBox, Select, fixed

def dict2tuple(d):
    return list(*zip(*zip(d.items())))

def inversedict(d):
    return {v:k for k, v in d.items()}

playerid_mapping = data['playerid_enc'].to_dict()
playerid_inverse = inversedict(playerid_mapping)

player_select = SelectMultiple(
    options=dict2tuple(playerid_mapping),
    value=(0,)
)

metric_select = Select(
    options=['p', 'pps'],
)

def scrub_posterior(player_encs: list, metric: str, trace):
    fig, ax = plt.subplots()
    for enc in player_encs:
        x, y = ECDF(trace[metric][:, enc])
        ax.plot(x, y, label=playerid_inverse[enc])
    ax.legend()
    ax.set_xlabel("p per salary")
    ax.set_ylabel('cumulative probability')
        
    sns.despine()
    
    
interact(scrub_posterior, player_encs=player_select, metric=metric_select, trace=fixed(trace_batting));

#### Hands-on: Finding our ideal player

Using the visualization, find a bunch of players whom you think would be great to scout, based on their PPS.

### Hierarchical model

**Discussion:** Is it reasonable that some players have a `p` (probability of hitting) value that ranges from as small as 0 to as high as 1? Under what assumption would this be reasonable? Under what assumption would this be *unreasonable*?

#### Code-along: Model Definition

To build a hierarchical model, we essentially need to define a prior on the parameters of the Beta distribution. These new priors now become "hyperpriors". How do we pick hyperpriors then? Some qualitative rules of thumb:

1. Pick a distribution with the correct **support**. For example, if a parameter can/should never take negative values, then pick a distribution with only positive support.
1. The distribution should be **weakly informative**. For example, rather than a uniform distribution hyperprior, one might choose a wide distribution (relative to the relevant scale of the parameters) that spans 1 order of magnitude difference, but still has uneven allocation of credibility. 

Use the PyMC3 distribution library to your advantage! There are pictures that show you the shapes of the distributions, which can be helpful in narrowing your choices.

with pm.Model() as hierarchical_baseball:
    a_prior = pm.Exponential("a_prior", lam=1/29.)
    b_prior = pm.Exponential("b_prior", lam=1/29.)
    
    p = pm.Beta("p", alpha=a_prior, beta=b_prior, shape=(len(data)))
    
    like = pm.Binomial("like", p=p, n=data['ab'], observed = data['h'])
    
    pps = pm.Deterministic("pps", p / data['salary'])

#### Sample from Posterior

with hierarchical_baseball:
    trace_hierarchical = pm.sample(2000)

While sampling happens, do you have any questions? Please feel free to ask.

#### Hands-on: Model Interpretation

interact(scrub_posterior, player_encs=player_select, metric=metric_select, trace=fixed(trace_hierarchical));

**Discussion:** 

1. What do you notice about the posterior distributions on `p`?
1. What do you notice about the posterior distributions on `pps`?

## Model Building Workflow

**Things we did:**

1. Get data into correct shape.
1. Build naive model, followed by more complex model.
1. Build a simple loss function into our model.


**Things we didn't do:**

1. Posterior predictive checks.

**Bonus things:**

1. Build custom posterior explorers. You can use your favourite framework of choice.





### Hands-on: Parameter estimation II -- the mean of a population

In this exercise, you'll calculate the  posterior mean beak depth of Galapagos finches in a given species. First you'll load the data and subset wrt species:

# Import and view head of data
df_12 = pd.read_csv('../data/finch_beaks_2012.csv')
df_fortis = df_12.loc[df_12['species'] == 'fortis']
df_scandens = df_12.loc[df_12['species'] == 'scandens']

To specify the full probabilty model, you need
- a likelihood function for the data &
- priors for all unknowns.

What is the likelihood here? Let's plot the measurements below and see that they look approximately Gaussian/normal so you'll use a normal likelihood $y_i\sim \mathcal{N}(\mu, \sigma^2)$. The unknowns here are the mean $\mu$ and standard deviation $\sigma$ and we'll use weakly informative priors on both
- a normal prior for $\mu$ with mean $10$ and standard deviation $5$;
- a uniform prior for $\sigma$ bounded between $0$ and $10$.

We can discuss biological reasons for these priors also but you can also test that the posteriors are relatively robust to the choice of prior here due to the amount of data.

sns.distplot(df_fortis['blength']);

with pm.Model() as model:
    # Prior for mean & standard deviation
    μ_1 = pm.Normal('μ_1', mu=10, sd=5)
    σ_1 = pm.Lognormal('σ_1', 0, 10)
    # Gaussian Likelihood
    y_1 = pm.Normal('y_1', mu=μ_1, sd=σ_1, observed=df_fortis['blength'])

# bust it out & sample
with model:
    samples = pm.sample(2000, njobs=1)

az.plot_posterior(samples);

## 2. Bayesian Hypothesis testing

### Bayesian Hypothesis testing I: A/B tests on click through rates

Assume we have a website and want to redesign the layout (*A*) and test whether the new layout (*B*) results in a higher click through rate. When people come to our website we randomly show them layout *A* or *B* and see how many people click through for each. First let's generate the data we need:

# click-through rates
p_a = 0.15
p_b = 0.20
N = 1000
n_successes_a = np.sum(np.random.uniform(size=N) <= p_a)
n_successes_b = np.sum(np.random.uniform(size=N) <= p_b)

Once again, we need to specify our models for $p_a$ and $p_b$. Each will be the same as the CTR example above
- Binomial likelihoods
- uniform priors on $p_a$ and $_p$.

We also want to calculate the posterior of the difference $p_a-p_b$ and we do so using `pm.Deterministic()`, which specifies a deterministic random variable, i.e., one that is completely determined by the values it references, in the case $p_a$ & $p_b$.

We'll now build the model:

with pm.Model() as Model:
    # Prior on p
    prob_a = pm.Uniform('p_a')
    prob_b = pm.Uniform('p_b')
    # Binomial Likelihood
    y_a = pm.Binomial('y_a', n=N, p=prob_a, observed=n_successes_a)
    y_b = pm.Binomial('y_b', n=N, p=prob_b, observed=n_successes_b)
    diff_clicks = pm.Deterministic('diff_clicks', prob_a-prob_b)

Sample from the posterior and plot them:

with Model:
    samples = pm.sample(2000, njobs=1)
az.plot_posterior(samples, kind='hist');

### Hands-on: Bayesian Hypothesis testing II -- beak lengths difference between species

**Task**: Determine whether the mean beak length of the Galapogas finches differs between species. For the mean of each species, use the same model as in previous hand-on section:

- Gaussian likelihood;
- Normal prior for the means;
- Uniform prior for the variances.

Also calculate the difference between the means and, for bonus points, the _effect size_, which is the difference between the means divided by the pooled standard deviations = $\sqrt{(\sigma_1^2+\sigma_2^2)/2}$. Hugo will talk through the importance of the _effect size_.

Don't forget to sample from the posteriors and plot them!

with pm.Model() as model:
    # Priors for means and variances
    μ_1 = pm.Normal('μ_1', mu=10, sd=5)
    σ_1 = pm.Uniform('σ_1', 0, 10)
    μ_2 = pm.Normal('μ_2', mu=10, sd=5)
    σ_2 = pm.Uniform('σ_2', 0, 10)
    # Gaussian Likelihoods
    y_1 = pm.Normal('y_1', mu=μ_1, sd=σ_1, observed=df_fortis['blength'])
    y_2 = pm.Normal('y_2', mu=μ_2, sd=σ_2, observed=df_scandens['blength'])
    # Calculate the effect size and its uncertainty.
    diff_means = pm.Deterministic('diff_means', μ_1 - μ_2)
    pooled_sd = pm.Deterministic('pooled_sd', 
                                 np.sqrt(np.power(σ_1, 2) + 
                                         np.power(σ_2, 2)) / 2)
    effect_size = pm.Deterministic('effect_size', 
                                   diff_means / pooled_sd)

# bust it out & sample
with model:
    samples = pm.sample(2000, njobs=1)
az.plot_posterior(samples, var_names=['μ_1', 'μ_2', 'diff_means', 'effect_size'], kind='hist');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08600817945222457, 'fraction_numerical': 0.011500805552113025, 'mean_word_length': 4.380184024536605, 'pattern_counts': {'":': 0, '<': 99, '<?xml version=': 0, '>': 95, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '621463', 'license': 'MIT', 'n_tokens_mistral': 12286, 'n_tokens_neox': 11472, 'n_words': 4956, 'path': 'notebooks/03-instructor-bayesian-estimation-win.ipynb', 'repo_name': 'masterfulEJ/bayesian-stats-modelling-tutorial', 'size': '40357'}