<jupyter_start><jupyter_text>IntroductionBaseball players have many metrics measured for them. Let's say we are on a baseball team, and would like to quantify player performance, one metric being their batting average (defined by how many times a batter hit a pitched ball, divided by the number of times they were up for batting ("at bat")). How would you go about this task?**Discuss**: 1. What data would we need?1. What metric would you rank by? 1. Would your metric be reasonable for rookie players?Let's read in the data.<jupyter_code>df = load_baseball()<jupyter_output><empty_output><jupyter_text>Betas, Bernoullis and Binomials: A brief introduction- Bernoulli distribution: a probability distribution modelling one "coin flip"-like trial. Parameterized by a single parameter `p`, which indicates probability of "success".- Binomial distribution: a probability distribution modelling the number of successes in `n` trials. Parameterized by both `n` and `p`.- Beta distributions: a probability distribution bounded over the interval $(0, 1)$. Models distribution of probability values, usually the `p` in a Bernoulli or Binomial. Parameterized by $\alpha$ and $\beta$, which can be thought of as "number of successes" and "number of failures" respectively.Every distribution has its "story". If you're curious, check out [<NAME>' probability stories][probstory] page.[probstory]: http://bebi103.caltech.edu.s3-website-us-east-1.amazonaws.com/2017/tutorials/t3b_probability_stories.htmlBeta-distribution Focus on betaLet's say we wanted to model a probability distribution centered approximately on 0.2. Depending on our parameterization of the Beta distribution, we can express different levels of confidence (as measured by the spread of the distribution) as to how sure we are a distribution takes on that value.The code + chart below should illustrate this clearly.<jupyter_code># To correctly express the beta distribution, we subtract 1.
low_prior = np.random.beta(a=1, b=7, size=100000)
med_prior = np.random.beta(a=19, b=79, size=100000)
high_prior = np.random.beta(a=199, b=799, size=100000)

fig = plt.figure()
ax = fig.add_subplot(111)

labels = ['low', 'med', 'high']
priors = [low_prior, med_prior, high_prior]

for label, prior in zip(labels, priors):
    x, y = ECDF(prior)
    ax.plot(x, y, label=label)
ax.legend()
ax.set_xlim(0, 1)
ax.set_xlabel('probability parameter')
ax.set_ylabel('cumulative fraction')
despine(ax)<jupyter_output><empty_output><jupyter_text>Notice how increasing the scale of parameterization (1x, 10x, 100x) decreased the variance. This leads to another way of looking at the Beta distribution:- $\alpha$ is the probability of success $p$ multiplied by a scale factor $\kappa$.- $\beta$ is $1-p$ multiplied by the same scale factor $\kappa$.We will see this in action below! **Exercise:** Write a naive estimation model for the players above.Hint, a possible model you could specify is as follows:![](../images/baseball-model.jpg)<jupyter_code>with pm.Model() as baseline_model:
    # Define a prior distribution for the probability parameter `p`, called
    # `thetas`, which is drawn from a Beta(0, 0) distribution, and is of
    # the same length as the 
    thetas = 
    # Define a likelihood function, which is the Binomial distribution,
    # in which `n` is the number of observed "at bats", p is `theta`, and 
    # the data we observe is the number of hits for that player.
    like =<jupyter_output><empty_output><jupyter_text>Now, sample from the posterior<jupyter_code># Your code below<jupyter_output><empty_output><jupyter_text>Plot the traceplots to check for convergence.<jupyter_code># Your code below.<jupyter_output><empty_output><jupyter_text>Visualize the posterior distribution using forest plots.<jupyter_code>ylabels = ylabels = "AB: " + df['AB'].astype(str) + ', H: ' + df['H'].astype('str')
# Your code below<jupyter_output><empty_output><jupyter_text>**Discuss:** Are the estimates reasonable, particularly for players that have had only one at bat (AB)?  Hierarchical Modelling**Discuss:** - How do we deal with the fact that some players have only had 1 at bat (AB = 1), and zero hits (H = 0)? - Would it be reasonable, fair, and in line with prior knowledge that the player's true batting average was zero?  AssumptionsIf it is qualitatively (and maybe also quantitatively) justifiable, we can impose the following assumption on the modelling process.We can assume that there is some underlying distribution for a player's batting average. We can then assume that each new player draws his own batting average from the underlying distribution. This encodes our intuition that human performance falls on some continuum that has a region of "concentrated" average. (Think "normal-like" distribution.)The modelling choices for the hierarchical model below are detailed in the [PyMC3 docs](https://docs.pymc.io/notebooks/hierarchical_partial_pooling.html), but are copied below for convenience. > We will assume that there exists a hidden factor (phi) related to the expected performance for all players (not limited to our 18). Since the population mean is an unknown value between 0 and 1, it must be bounded from below and above. Also, we assume that nothing is known about global average. Hence, a natural choice for a prior distribution is the uniform distribution.> > Next, we introduce a hyperparameter kappa to account for the variance in the population batting averages, for which we will use a bounded Pareto distribution. This will ensure that the estimated value falls within reasonable bounds. These hyperparameters will be, in turn, used to parameterize a beta distribution, which is ideal for modeling quantities on the unit interval. The beta distribution is typically parameterized via a scale and shape parameter, it may also be parametrized in terms of its mean $\mu \in [0, 1]$ and sample size (a proxy for variance) $\nu = \alpha + \beta (\nu > 0)$>> The final step is to specify a sampling distribution for the data (hit or miss) for every player, using a Binomial distribution. This is where the data are brought to bear on the model.> > We could use `pm.Pareto('kappa', m=1.5)`, to define our prior on kappa, but the Pareto distribution has very long tails. Exploring these properly is difficult for the sampler, so we use an equivalent but faster parametrization using the exponential distribution. We use the fact that the log of a Pareto distributed random variable follows an exponential distribution.The most important thing to note is how the ${\theta}$ parameters are drawn from a parental distribution that is not broadcast across the samples.**Notes:** - If you recall the parameterization of the beta distribution described above, it is exactly $p$ and $1-p$ multiplied by a scalar $\kappa$! In the model spec below, we are treating phi and kappa as being parameters to estimate.- The model below expresses the intuition that the scale parameter is Pareto distributed - the intuition being that approx. 80% of the distribution is going to be small, but 20% of the time it'll be large.In visual form, this is the model below.![](../images/baseball-hierarchical-model.jpg)<jupyter_code>with pm.Model() as baseball_model:
    
    phi = pm.Uniform('phi', lower=0.0, upper=1.0)
    kappa_log = pm.Exponential('kappa_log', lam=1.5)
    kappa = pm.Deterministic('kappa', tt.exp(kappa_log))
    
    thetas = pm.Beta('thetas', alpha=phi*kappa, beta=(1.0-phi)*kappa, shape=len(df))
    like = pm.Binomial('like', n=df['AB'], p=thetas, observed=df['H'])<jupyter_output><empty_output><jupyter_text>Sample from the posterior.<jupyter_code># Your code below<jupyter_output><empty_output><jupyter_text>Visualize the trace plots to check for convergence.<jupyter_code># Your code below<jupyter_output><empty_output><jupyter_text>Visualize the posterior distributions using forestplots.<jupyter_code>ylabels = "AB: " + df['AB'].astype(str) + ', H: ' + df['H'].astype('str')
# Your code below.<jupyter_output><empty_output><jupyter_text>With a hierarchical model, we make the assumption that our observations (or treatments that group our observations) are somehow related. Under this assumption, when we have a new sample for which we have very few observations, we are able to borrow power from the population to make inferences about the new sample. Depending on the scenario, this assumption can either be reasonable, thereby not necessitating much debate, or be considered a "strong assumption", thereby requiring strong justification.   Shrinkage"Shrinkage" is a term used to describe how hierarchical model estimation will usually result in parameter estimates that are "shrunk" away from their maximum likelihood estimators (i.e. the naive estimate from the data) towards the global mean. Shrinkage in and of itself is not necessarily a good or bad thing. However, because hierarchical models can sometimes be tricky to get right, we can use a shrinkage plot as a visual diagnostic for whether we have implemented the model correctly.<jupyter_code># MLE per player
no_pool = df['batting_avg'].values
# Hierarchical model
partial_pool = trace['thetas'].mean(axis=0)
# MLE over all players
complete_pool = np.array([df['batting_avg'].mean()] * len(df))

fig = plt.figure()
ax = fig.add_subplot(111)

for r in np.vstack([no_pool, partial_pool, complete_pool]).T:
    ax.plot(r)
    
ax.set_xticks([0, 1, 2])
ax.set_xticklabels(['no pooling', 'partial pooling', 'complete pooling'])
ax.set_ylim(0, 1)
despine(ax)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08418073173288605, 'fraction_numerical': 0.010168780794632562, 'mean_word_length': 5.2109375, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 54, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '218028', 'license': 'MIT', 'n_tokens_mistral': 2713, 'n_tokens_neox': 2537, 'n_words': 1215, 'path': 'notebooks/06-student-hierarchical-baseball.ipynb', 'repo_name': 'IsabelAverill/Scipy-2018---Bayesian', 'size': '9544'}