# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Think Bayes
#
# Second Edition
#
# Copyright 2020 <NAME>
#
# License: [Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)](https://creativecommons.org/licenses/by-nc-sa/4.0/)

# +
# If we're running on Colab, install empiricaldist
# https://pypi.org/project/empiricaldist/

import sys
IN_COLAB = 'google.colab' in sys.modules

if IN_COLAB:
    # !pip install empiricaldist

# +
# Get utils.py and create directories

import os

if not os.path.exists('utils.py'):
    # !wget https://github.com/AllenDowney/ThinkBayes2/raw/master/code/soln/utils.py
        
if not os.path.exists('figs'):
    # !mkdir figs
        
if not os.path.exists('tables'):
    # !mkdir tables

# +
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from empiricaldist import Pmf, Cdf
from utils import decorate, savefig, write_table
# -

# ## Bayesian inference
#
# Whenever people compare the Bayesian inference with conventional approaches, one of the questions that comes up most often is something like, "What about p-values?"
# And one of the most common examples is the comparison of two groups to see if there is a difference in their means.
#
# In classical statistical inference, the usual tool for this scenario is a [Student's *t*-test](https://en.wikipedia.org/wiki/Student%27s_t-test), and the result is a [p-value](https://en.wikipedia.org/wiki/P-value).
# This process is an example of [null hypothesis significance testing](https://en.wikipedia.org/wiki/Statistical_hypothesis_testing).
#
# A Bayesian alternative is to compute the posterior distribution of the difference between the groups.
# Then we can use that distribution to answer whatever questions we are interested in, including the most likely size of the difference, a credible interval that's likely to contain the true difference, the probability of superiority, or the probability that the difference exceeds some threshold.
#
# To demonstrate this process, I'll solve a standard problem from a statistical textbook, comparing the effect of an educational "treatment" compared to a control.

# ## Improving Reading Ability
#
# We'll use data from a [Ph.D. dissertation in educational psychology](https://docs.lib.purdue.edu/dissertations/AAI8807671/) written in 1987, which was used as an example in a [statistics textbook](https://books.google.com/books/about/Introduction_to_the_practice_of_statisti.html?id=pGBNhajABlUC) from 1989 and published on [DASL](https://web.archive.org/web/20000603124754/http://lib.stat.cmu.edu/DASL/Datafiles/DRPScores.html), a web page that collects data stories.  
#
# Here's the description from DASL:
#
# > An educator conducted an experiment to test whether new directed reading activities in the classroom will help elementary school pupils improve some aspects of their reading ability. She arranged for a third grade class of 21 students to follow these activities for an 8-week period. A control classroom of 23 third graders followed the same curriculum without the activities. At the end of the 8 weeks, all students took a Degree of Reading Power (DRP) test, which measures the aspects of reading ability that the treatment is designed to improve.
#
# The [dataset is available here](https://web.archive.org/web/20000603124754/http://lib.stat.cmu.edu/DASL/Datafiles/DRPScores.html); I have put it in a CSV file, which we can download like this:

# +
import os

if not os.path.exists('drp_scores.csv'):
    # !wget https://github.com/AllenDowney/ThinkBayes2/raw/master/data/drp_scores.csv
# -

# I'll use Pandas to load the data into a `DataFrame`.

# +
import pandas as pd

df = pd.read_csv('drp_scores.csv', skiprows=21, delimiter='\t')
df.head()
# -

# I'll use `groupby` to separate the data for the `Treated` and `Control` groups:

# +
grouped = df.groupby('Treatment')
responses = {}

for name, group in grouped:
    responses[name] = group['Response']
# -

# Here are CDFs of the scores for the two groups and summary statistics.

# +
for name, response in responses.items():
    print(name, len(response), response.mean(), response.std())
    cdf = Cdf.from_seq(response)
    cdf.plot(label=name)
    
decorate(xlabel='Score', ylabel='CDF')
savefig('fig11-01')
# -

# The distribution of scores is not exactly normal for either group, but it is close enough that the normal model is a reasonable choice.
#
# So I'll assume that in the entire population of students (not just the ones in the experiment), the distribution of scores is well modeled by a normal distribution with unknown mean and standard deviation.
# I'll use `mu` and `sigma` to denote these unknown population parameters.
#
# And we'll do a Bayesian update to estimate what they are.

# ## Estimating parameters
#
# As always, we need a prior distribution for the parameters.  
# Since there are two parameters, it will be a joint distribution.  
# I'll construct it by choosing marginal distributions for each parameter and computing their outer product.
#
# As a simple starting place, I'll assume that the prior distributions for `mu` and `sigma` are uniform.

mus = np.linspace(20, 80, 101)
prior_mu = Pmf(1, mus, name='mean')

sigmas = np.linspace(5, 30, 101)
prior_sigma = Pmf(1, sigmas, name='std')

# +
from utils import outer_product

prior = outer_product(prior_mu, prior_sigma)
# -

data = responses['Control']
data.shape

# Now, we would like to know the probability of each score in the dataset for each hypothetical pair of values, `mu` and `sigma`.
# I'll do that by making a 3-dimensional grid with values of `sigma` on the first axis, values of `mu` on the second axis, and the scores from the dataset on the third axis.

sigmas, mus, data_mesh = np.meshgrid(prior.columns, prior.index, data)
sigmas.shape

# Now we can use `norm.pdf` to compute the probability density of each score for each hypothetical pair of parameters.

# +
from scipy.stats import norm

densities = norm.pdf(data_mesh, sigmas, mus)
densities.shape
# -

# The result is a 3-D array.  To compute likelihoods, I'll compute the product of these densities along the third axis, that is `axis=2`:

likelihood = densities.prod(axis=2)
likelihood.shape

# The result is a 2-D array that contains the likelihood of the entire dataset for each hypothetical pair of parameters.
#
# We can use this array as part of a Bayesian update, as in this function:

# +
from utils import normalize

def update_norm(prior, data):
    """Update the prior based on data.
    
    prior: joint distribution of mu and sigma
    data: sequence of observations
    """
    X, Y, Z = np.meshgrid(prior.columns, prior.index, data)
    likelihood = norm.pdf(Z, Y, X).prod(axis=2)

    posterior = prior * likelihood
    normalize(posterior)

    return posterior


# -

# Here are the updates for the control and treatment groups:

data = responses['Control']
posterior_control = update_norm(prior, data)

data = responses['Treated']
posterior_treated = update_norm(prior, data)


# And here's what they look like:

def plot_contour(joint, **options):
    """Plot a joint distribution.
    
    joint: DataFrame representing a joint PMF
    """
    cs = plt.contour(joint.columns, joint.index, joint, **options)
    decorate(xlabel='Standard deviation', ylabel='Mean')
    return cs


# +
plot_contour(posterior_control, cmap='Blues')
plt.text(18, 49.5, 'Control', color='C0')

cs = plot_contour(posterior_treated, cmap='Oranges')
plt.text(12, 57, 'Treated', color='C1')

savefig('fig11-02')
# -

# Along the vertical axis, it looks like the mean score for the treated group is higher.
# Along the horizontal axis, it looks like the standard deviation for the control group is higher.
#
# If we think the treatment causes these differences, the data suggest that the treatment increases the mean score and decreases their spread.
# We can see these differences more clearly by looking at the marginal distributions for `mu` and `sigma`.

# ## Posterior marginal distributions
#
# I'll use `marginal`, which we saw in Chapter 9, to extract the posterior marginal distributions for the population means.

# +
from utils import marginal

pmf_mean_control = marginal(posterior_control, 1)
pmf_mean_treated = marginal(posterior_treated, 1)
# -

# Here's what they look like:

# +
pmf_mean_control.plot(label='Control')
pmf_mean_treated.plot(label='Treated')

decorate(xlabel='Population mean', 
         ylabel='PDF', 
         title='Posterior distributions of mu')

savefig('fig11-03')
# -

# It looks like we are pretty sure that the population mean in the treated group is higher.
# We can use `prob_gt` to compute the probability of superiority:

Pmf.prob_gt(pmf_mean_treated, pmf_mean_control)

# There is a 98% chance that the mean in the treated group is higher.
#
# We can use `sub_dist` to compute the distribution of the difference.

diff = Pmf.sub_dist(pmf_mean_treated, pmf_mean_control)

# Two things to be careful about when you use methods like `sub_dist`: 
#
# The first is that the result usually contains more elements than the original `Pmf`.  
# In this example, the original distributions have the same quantities, so the size increase is moderate.

len(pmf_mean_treated), len(pmf_mean_control), len(diff)

# In the worst case, the size of the result can be the product of the sizes of the originals.
#
# The other thing to be aware of is that plotting a `Pmf` does not always work well.
# In this example, if we plot the distribution of differences, the result is pretty noisy.

# +
diff.plot()

decorate(xlabel='Difference in population means', 
         ylabel='PDF', 
         title='Posterior distribution of difference in mu')
# -

# There are two ways to work around that limitation.  One is to plot the CDF, which smooths out the noise:

# +
diff.make_cdf().plot()

decorate(xlabel='Difference in population means', 
         ylabel='CDF', 
         title='Posterior distribution of difference in mu')
# -

# The other option is to use kernel density estimation (KDE) to make a smooth approximation of the PDF on an equally-spaced grid.

# +
from scipy.stats import gaussian_kde

def make_kde(pmf, n=101):
    """Make a kernel density estimate for a PMF.
    
    pmf: Pmf object
    n: number of points
    
    returns: Pmf object
    """
    kde = gaussian_kde(pmf.qs, weights=pmf.ps)
    qs = np.linspace(pmf.qs.min(), pmf.qs.max(), n)
    ps = kde.evaluate(qs)
    pmf = Pmf(ps, qs)
    pmf.normalize()
    return pmf


# -

# Here's what it looks like.

# +
make_kde(diff).plot()

decorate(xlabel='Difference in population means', 
         ylabel='PDF', 
         title='Posterior distribution of difference in mu')

savefig('fig11-04')
# -

# The mean is almost 10 points, which is substantial.

diff.mean()

# Finally, we can use `credible_interval` to compute a 90% credible interval.

diff.credible_interval(0.9)

# Based on the data, we are pretty sure the treatment improves test scores by 2 to 17 points.

# ## Using summary statistics
#
# In this example the dataset is not very big, so it doesn't take too long to compute the probability of every score under every hypothesis.
# But the result is a 3-D array; for larger datasets, it might be too big to compute practically.
#
# Also, with larger datasets the likelihoods get very small, sometimes so small that we can't compute them with normal floating-point arithmetic.
# That's because we are computing the probability of a particular dataset; the number of possible datasets is astronomically big, so the probability of any of them is very small.
#
# An alternative is to compute a summary of the dataset and compute the likelihood of the summary.
# For example, if we compute the sample mean of the data and the sample standard deviation, we could compute the likelihood of those summary statistics under each hypothesis.
#
# As an example, suppose we know that the population mean is 40 and the standard deviation is 17.
# We can make a `norm` object that represents a normal distribution with these parameters:

mu = 40
sigma = 17
dist = norm(mu, sigma)

# Now suppose we draw 1000 samples from this distribution with sample size `n=20`.
# I'll use `rvs`, which generates a random sample, to simulate this experiment.

n = 20
samples = dist.rvs((1000, n))
samples.shape

# The result is an array with 1000 rows, each containing a sample with 20 columns.
#
# If we compute the mean of each row, the result is an array that contains 1000 sample means; that is, each value is the mean of a sample with `n=20`.

sample_means = samples.mean(axis=1)
sample_means.shape

# Now, we would like to know what the distribution of these sample means is.
# Using the properties of the normal distribution, [we can show](https://en.wikipedia.org/wiki/Sum_of_normally_distributed_random_variables) that their distribution is normal with mean $\mu$ and standard deviation $\sigma/\sqrt{n}$:

dist_m = norm(mu, sigma/np.sqrt(n))

# `dist_m` represents the "sampling distribution of the mean".  
# We can use it to make a `Cdf` that approximates the CDF of the sampling distribution

low, high = sample_means.min(), sample_means.max()
qs = np.linspace(low, high, 101)
ps = dist_m.cdf(qs)
cdf_m = Cdf(ps, qs)

# The following figure shows this theoretical distribution along with the empirical distribution of the values in `sample_means`.

# +
cdf_m.plot(label='Sampling distribution of the mean')
Cdf.from_seq(sample_means).plot(label='Distribution of sample means')

decorate(xlabel='Mean score',
         ylabel='CDF')
# -

# The random sample means follow the theoretical distribution closely, as expected.
#
# We can also compute standard deviations for each row in `samples`.

sample_stds = samples.std(axis=1)
sample_stds.shape

# The result is an array of sample standard deviations.
# We might wonder what the distribution of these values is.
# The [derivation is not as easy](https://en.wikipedia.org/wiki/Normal_distribution#Sample_variance), but if we transform the sample standard deviations like this:
#
# $t = n s^2 / \sigma^2$
#
# where $n$ is the sample size, $s$ is the sample standard deviation, and $\sigma$ is the population standard deviation, the transformed values follow a [chi-square distribution](https://en.wikipedia.org/wiki/Chi-square_distribution) with $n-1$ degrees of freedom.
#
# Here are the transformed values.

transformed = n * sample_stds**2 / sigma**2

# And I'll create a `chi2` object that represents a chi-square distribution.

# +
from scipy.stats import chi2

dist_s = chi2(n-1)
# -

# And a `Cdf` that approximates the CDF of the chi-square distribution.

low, high = transformed.min(), transformed.max()
qs = np.linspace(low, high, 101)
ps = dist_s.cdf(qs)
cdf_s = Cdf(ps, qs)

# Now we can compare the theoretical distribution to the empirical distribution of the sample standard deviations.

# +
cdf_s.plot(label='Sampling distribution of s')
Cdf.from_seq(transformed).plot(label='Distribution of sample s')

decorate(xlabel='Standard deviation of scores (transformed)',
         ylabel='CDF')
# -

# The distribution of transformed sample standard deviations agrees with the theoretical distribution.
#
# I think it is useful to check theoretical results like this, for a few reasons:
#
# * It confirms that my understanding of the theory is correct,
#
# * It confirms that the conditions where I am applying the theory are conditions where the theory holds, 
#
# * It confirms that the implementation details are correct.  For many distributions, there is more than one way to specify the parameters.  If you use the wrong specification, this kind of testing will help you catch the error.

# Before we move on, I'll mention one other theoretical result we will use: [Basu's theorem](https://en.wikipedia.org/wiki/Basu%27s_theorem), which states that the sample mean and sample standard deviation are independent.
#
# To see whether that's true, we can make a contour plot for the joint distribution of sample mean and sample standard deviation.

# +
import seaborn as sns

sns.kdeplot(sample_stds, sample_means)

decorate(xlabel='Sample standard deviation', ylabel='Sample mean')
# -

# It looks like the axes of the ellipses are aligned with the axes, which indicates that the variables are independent.

# ## Update with summary statistics
#
# Now we're ready to do an update.
# I'll compute summary statistics for the two groups.

# +
summary = {}

for name, response in responses.items():
    summary[name] = len(response), response.mean(), response.std()
    
summary
# -

# The result is a dictionary that maps from group name to a tuple that contains the sample size, `n`, the sample mean, `m`, and the sample standard deviation `s`, for each group.
#
# I'll demonstrate the update with the summary statistics from the control group.

n, m, s = summary['Control']

# I'll make a mesh with hypothetical values of `mu` on the vertical axis and values of `sigma` on the horizontal axis.

sigmas, mus = np.meshgrid(prior.columns, prior.index)
sigmas.shape

# Now we can compute the likelihood of seeing the sample mean `m` for each pair of parameters.

like1 = norm.pdf(m, mus, sigmas/np.sqrt(n))

# And use it to update the prior.

posterior1 = prior * like1
normalize(posterior1)
plot_contour(posterior1);

# The joint posterior distribution is narrow along the vertical axis and wide along the horizontal axis, which indicates that the sample mean provides a lot of information about the population mean, but by itself it doesn't tell us much about the population standard deviation.
#
# Let's see what happens if we update with only the standard deviation and not the mean.
# Here's the likelihood:

like2 = chi2.pdf(n * s**2 / sigmas**2, n-1)

# And here's the update:

posterior2 = prior * like2
normalize(posterior2)
plot_contour(posterior2);


# The posterior joint distribution is narrow along the horizontal axis and very wide along the vertical axis.
#
# So the sample standard deviation provides information about the population standard deviation, but no information at all about the mean.
#
# The following function does both updates, using the sample mean and standard deviation.

def update_norm_summary(prior, data):
    """Update a normal distribution based on summary statistics.
    
    prior: DataFrame, joint prior distribution
    data: tuple of sample size, sample mean, sample std
    """
    n, m, s = data
    sigmas, mus = np.meshgrid(prior.columns, prior.index)
    
    like1 = norm.pdf(m, mus, sigmas/np.sqrt(n))
    like2 = chi2.pdf(n * s**2 / sigmas**2, n-1)
    
    posterior = prior * like1 * like2
    normalize(posterior)
    
    return posterior


# Here are the updates for the two groups.

data = summary['Control']
posterior_control2 = update_norm_summary(prior, data)

data = summary['Treated']
posterior_treated2 = update_norm_summary(prior, data)

# And here are the results.

# +
plot_contour(posterior_control2, cmap='Blues')
plt.text(18, 49.5, 'Control', color='C0')

cs = plot_contour(posterior_treated2, cmap='Oranges')
plt.text(12, 57, 'Treated', color='C1');
# -

# Visually, these posterior joint distributions are similar to the ones we computed using the entire datasets, not just the summary statistics.
# But they are not exactly the same, as we can see by comparing the marginal distributions.

# ## Comparing marginals
#
# Again, let's extract the marginal posterior distributions.

# +
from utils import marginal

pmf_mean_control2 = marginal(posterior_control2, 1)
pmf_mean_treated2 = marginal(posterior_treated2, 1)
# -

# And compare them to results we got using the entire dataset.

# +
pmf_mean_control.plot(color='gray')
pmf_mean_control2.plot(label='Control')
pmf_mean_treated.plot(color='gray')
pmf_mean_treated2.plot(label='Treated')

decorate(xlabel='Population mean', 
         ylabel='PDF', 
         title='Posterior distributions of mu')

savefig('fig11-05')
# -

# For both groups, the distribution of `mu` is a little wider when we use only the summary statistics; that is, we are a little less certain about the values of the means.
#
# If we compute the posterior distribution of the difference in means, like this:

diff2 = Pmf.sub_dist(pmf_mean_treated2, pmf_mean_control2)
diff2.make_cdf().plot()

# The mean difference is nearly the same.

diff.mean(), diff2.mean()

# But the credible interval is a bit wider.

diff.credible_interval(0.9), diff2.credible_interval(0.9)

# That's because the update we did is based on the implicit assumption that the distribution of the data is actually normal.
# But it's not; as a result, when we replace the dataset with the summary statistics, we lose some information about the true distribution of the data.
# With less information, we are less certain about the parameters.

# ## Summary
#
# In this chapter we used a joint distribution to represent prior probabilities for the parameters of a normal distribution, `mu` and `sigma`.
#
# And we updated that distribution two ways: first using the entire dataset and the normal PDF; then using summary statistics, the normal PDF, and the chi-square PDF.
#
# Using summary statistics is computationally more efficient, but it loses some information in the process.
#
# Normal distributions appear in many domains, as well as other distributions that are well approximated by normal distributions.  So the methods in this chapter are broadly applicable.  The exercises at the end of the chapter will give you a chance to apply them.

# ## Exercises

# **Exercise:** Looking again at the posterior joint distribution of `mu` and `sigma`, it seems like the standard deviation of the treated group might be lower; if so, that would suggest that the treatment is more effective for students with lower scores.
#
# But before we speculate too much, we should estimate the size of the difference and see whether it might actually be 0.
#
# As we did with the values of `mu` in the previous section, extract the posterior marginal distributions of `sigma` for the two groups.
# What is the probability that the standard deviation is higher in the control group?
#
# Compute the distribution of the difference in `sigma` between the two groups.  What is the mean of this difference?  What is the 90% credible interval?

# +
# Solution

pmf_std_control = marginal(posterior_control, 0)
pmf_std_treated = marginal(posterior_treated, 0)

# +
# Solution

pmf_std_control.plot(label='Control')
pmf_std_treated.plot(label='Treated')

decorate(xlabel='Population standard deviation', ylabel='PDF', 
         title='Posterior distributions of sigma')

# +
# Solution

Pmf.prob_gt(pmf_std_control, pmf_std_treated)

# +
# Solution

diff = Pmf.sub_dist(pmf_std_control, pmf_std_treated)

# +
# Solution

diff.mean()

# +
# Solution

diff.credible_interval(0.9)

# +
# Solution

make_kde(diff).plot()

decorate(xlabel='Difference in population standard deviation', 
         ylabel='PDF', 
         title='Posterior distributions of difference in sigma')
# -

# ## Exercise
#
# An "[effect size](http://en.wikipedia.org/wiki/Effect_size)" is a statistic intended to quantify the magnitude of a phenomenon.
# If the phenomenon is a difference in means between two groups, a common way to quantify it is Cohen's effect size, denoted $d$.
#
# If the parameters for Group 1 are $(\mu_1, \sigma_1)$, and the
# parameters for Group 2 are $(\mu_2, \sigma_2)$, Cohen's
# effect size is
#
# $ d = \frac{\mu_1 - \mu_2}{(\sigma_1 + \sigma_2)/2} $
#
# Use the joint posterior distributions for the two groups to compute the posterior distribution for Cohen's effect size.
#
# Hint: if enumerating all pairs from the two distributions takes too
# long, consider random sampling.  One way to do that is to "stack" the joint distribution and convert it to a `Pmf`.

pmf = Pmf(posterior_treated.stack())
pmf.qs

# The result is a `Pmf` where the elements of `qs` are tuples containing possible pairs of $\mu$ and $\sigma$.  The `ps` are the corresponding probabilities for each pair.
#
# `Pmf` provides `choice`, which we can use to draw a sample of parameters from the posterior distribution.

sample_treated = pmf.choice(1000)
param1 = sample_treated[0]
param1

# If you do the same for `posterior_control`, you can use the two samples to estimate the distribution of Cohen's effect size.  Then compute the mean and 90% credible interval. 

# +
# Solution

pmf = Pmf(posterior_control.stack())
sample_control = pmf.choice(1000)
param2 = sample_control[0]
param2


# +
# Solution

def cohen_effect(t):
    """Compute Cohen's effect size for difference in means.
    
    t: tuple of (mu1, sigma1), (mu2, sigma2)
    
    return: float
    """
    (mu1, sigma1), (mu2, sigma2) = t
    sigma = (sigma1 + sigma2) / 2
    return (mu1 - mu2) / sigma


# +
# Solution

cohen_effect((param1, param2))

# +
# Solution

a = np.transpose([sample_treated, sample_control])

ds = np.apply_along_axis(cohen_effect, 1, a)
ds.shape

# +
# Solution

cdf = Cdf.from_seq(ds)
cdf.plot()

decorate(xlabel='Cohen effect size',
         ylabel='CDF',
         title='Posterior distributions of effect size')

# +
# Solution

cdf.mean()

# +
# Solution

cdf.credible_interval(0.9)
# -

# **Exercise:** This exercise is inspired by [a question that appeared on Reddit](https://www.reddit.com/r/statistics/comments/hcvl2j/q_reverse_empirical_distribution_rule_question/).
#
# An instructor announces the results of an exam like this, "The average score on this exam was 81.  Out of 25 students, 5 got more than 90, and I am happy to report that no one failed (got less than 60)."
#
# Based on this information, what do you think the standard deviation of scores was?
#
# You can assume that the distribution of scores is approximately normal.  And let's assume that the sample mean, 81, is actually the population mean, so we only have to estimate `sigma`.
#
# Hint: To compute the probability of a score greater than 90, you can use `norm.sf`, which computes the survival function, also known as the complementary CDF, or `1 - cdf(x)`.

# +
# Solution

hypos = np.linspace(1, 41, 101)

# +
# Solution

# Here's are the probabilities of a score greater than 90
# for each hypothetical value of sigma.

from scipy.stats import norm

pgt90 = norm(81, hypos).sf(90)
pgt90.shape

# +
# Solution

# And here's the chance that 5 out of 25 people
# get a score greater than 90

from scipy.stats import binom

likelihood1 = binom(25, pgt90).pmf(5)
likelihood1.shape

# +
# Solution

prior = Pmf(1, hypos)
posterior = prior * likelihood1
posterior.normalize()

# +
# Solution

# Here's the posterior after the first update.

posterior.plot()
decorate(xlabel='Standard deviation',
         ylabel='PMF')

# +
# Solution

# Here's the probability of a score greater than 60

pgt60s = norm(81, hypos).sf(60)

# +
# Solution

# And here's the probability that all 25 students exceed 60

likelihood2 = pgt60s ** 25

# +
# Solution

plt.plot(hypos, likelihood2)
decorate(xlabel='Standard deviation',
         ylabel='Likelihood')

# +
# Solution

# Here's the posterior after both updates

prior = Pmf(1, hypos)
prior.normalize()
posterior2 = prior * likelihood1 * likelihood2
posterior2.normalize()

# +
# Solution

posterior.plot(label='Posterior 1')
posterior2.plot(label='Posterior 2')

decorate(xlabel='Standard deviation',
         ylabel='PMF')

# +
# Solution

posterior.mean(), posterior2.mean()

# +
# Solution

posterior2.credible_interval(0.9)


# -

# **Exercise:** I have a soft spot for crank science, so this exercise is about the [Variability Hypothesis](http://en.wikipedia.org/wiki/Variability_hypothesis), which
#
# > "originated in the early nineteenth century with <NAME>, who argued that males have a greater range of ability than females, especially in intelligence. In other words, he believed that most geniuses and most mentally retarded people are men. Because he considered males to be the ’superior animal,’ Meckel concluded that females’ lack of variation was a sign of inferiority."
#
# I particularly like that last part because I suspect that if it turned out that women were _more_ variable, Meckel would have taken that as a sign of inferiority, too. 
#
# Nevertheless, the Variability Hypothesis suggests an exercise we can use to practice the methods in this chapter.  Let's look at the distribution of heights for men and women in the U.S. and see who is more variable.
#
# I used 2018 data from the CDC’s [Behavioral Risk Factor Surveillance System](https://www.cdc.gov/brfss/annual_data/annual_2018.html) (BRFSS), which includes self-reported heights from 154407 men and 254722 women. 

# Here's what I found:
#
# * The average height for men is 178 cm; the average height for women is 163 cm. So men are taller on average; no surprise there.
#
# * For men the standard deviation is 8.27 cm; for women it is 7.75 cm. So in absolute terms, men's heights are more variable.
#
# But to compare variability between groups, it is more meaningful to use the [coefficient of variation](https://en.wikipedia.org/wiki/Coefficient_of_variation) (CV), which is the standard deviation divided by the mean. It is a dimensionless measure of variability relative to scale. 
#
# For men CV is 0.0465; for women it is 0.0475.
# The coefficient of variation is higher for women, so this dataset provides evidence against the Variability Hypothesis. But we can use Bayesian methods to make that conclusion more precise.

# Use these summary statistics to compute the posterior distribution of `mu` and `sigma` for the distributions of male and female height.
# Use `Pmf.div_dist` to compute posterior distributions of CV.
# Based on this dataset and the assumption that the distribution of height is normal, what is the probability that the coefficient of variation is higher for men?
# What is the most likely ratio of the CVs and what is the 90% credible interval for that ratio?
#
# Hint: Use different prior distributions for the two groups, and chose them so they cover all parameters with non-negligible probability.
#
# Also, you might find this function helpful:

def get_posterior_cv(joint):
    """Get the posterior distribution of CV.
    
    joint: joint distribution of mu and sigma
    
    returns: Pmf representing the smoothed posterior distribution
    """
    pmf_mean = marginal(joint, 1)
    pmf_std = marginal(joint, 0)
    pmf_cv = Pmf.div_dist(pmf_std, pmf_mean)
    return make_kde(pmf_cv)


# +
# Solution

n = 154407
mean = 178
std = 8.27

# +
# Solution

mus = np.linspace(mean-0.1, mean+0.1, 101)
prior_mu = Pmf(1, mus, name='mean')

sigmas = np.linspace(std-0.1, std+0.1, 101)
prior_sigma = Pmf(1, sigmas, name='std')

prior = outer_product(prior_mu, prior_sigma)

# +
# Solution

data = n, mean, std
posterior_male = update_norm_summary(prior, data)
plot_contour(posterior_male, cmap='Blues');

# +
# Solution

n = 254722
mean = 163
std = 7.75

# +
# Solution

mus = np.linspace(mean-0.1, mean+0.1, 101)
prior_mu = Pmf(1, mus, name='mean')

sigmas = np.linspace(std-0.1, std+0.1, 101)
prior_sigma = Pmf(1, sigmas, name='std')

prior = outer_product(prior_mu, prior_sigma)

# +
# Solution

data = n, mean, std
posterior_female = update_norm_summary(prior, data)
plot_contour(posterior_female, cmap='Oranges');

# +
# Solution

pmf_cv_male = get_posterior_cv(posterior_male)
make_kde(pmf_cv_male).plot()

pmf_cv_female = get_posterior_cv(posterior_female)
make_kde(pmf_cv_female).plot()

decorate(xlabel='Coefficient of variation',
         ylabel='PDF',
         title='Posterior distributions of CV')

# +
# Solution

ratio_cv = Pmf.div_dist(pmf_cv_female, pmf_cv_male)
ratio_cv.max_prob()

# +
# Solution

ratio_cv.credible_interval(0.9)
# -



StarcoderdataJupyterScriptsDedupFiltered