<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Poincaré's bread
Allen Downey
MIT License
Step1:  Poincaré's bread
Henri Poincaré was a French mathematician who taught at the Sorbonne around 1900.The following anecdote about him is probably fabricated, but it makes an interesting probability problem.
Supposedly, Poincaré suspected that his local bakery was selling loaves of bread that were lighter than the advertised weight of 1 kg, so every day for a year he bought a loaf of bread, brought it home and weighed it. At the end of the year, he plotted the distribution of his measurements and showed that it fit a normal distribution with mean 950 g and standard deviation 50 g. He brought this evidence to the bread police, who gave the baker a warning.
For the next year, Poincaré continued the practice of weighing his bread every day. At the end of the year, he found that the average weight was 1,000 g, just as it should be, but again he complained to the bread police, and this time they fined the baker.
Why? Because the shape of the distribution was asymmetric. Unlike the normal distribution, it was skewed to the right, which is consistent with the hypothesis that the baker was still making 950 g loaves, but deliberately giving Poincaré the heavier ones.
Write a program that simulates a baker who chooses n loaves from a distribution with mean 950 g and standard deviation 50 g, and gives the heaviest one to Poincaré. What value of n yields a distribution with mean 1000 g? What is the standard deviation?
Compare this distribution to a normal distribution with the same mean and the same standard deviation. Is the difference in the shape of the distribution big enough to convince the bread police?
First, I'll make a function that simulates choosing n loaves and giving Poincaré the heaviest.
Step2:  Here's a sample from the distribution of bread weight if the baker chooses 5 loaves.
Step3:  And here's how the observed mean depends on n.
Step4:  Looks like n=4 yields a measured mean close to 1000 g.
Here's a large sample of bread weights with n=4 so we can estimate the mean and standard deviation.
Step5:  Now let's compare the distribution of the sample to the normal distribution with the same mean and std.
Step6:  Visually they are almost indistinguishable, which makes it seem unlikely that Poincaré could really tell the difference, especially with only 365 observations.
But let's look more closely.  First, let's simulate one year.
Step7:  As a test statistic, we'll use the sample skew.
Step8:  The null hypothesis is that the data come from a normal distribution with the observed mean and standard deviation.
Step9:  Now we can generate fake data under the null hypothesis and compute the test statistic.
Step10:  And we can estimate the distribution of the test statistic under the null hypothesis
Step11:  Here's the 95th percentile of the test statistic under the null hypothesis.
Step12:  If Herni observes a skew above this threshold, he would conclude that it is statistically significant.
Let's see how often that would happen under the alternate hypothesis, that they baker is giving him the heaviest of 4 loaves.
Step13:  Here's what the sampling distribution looks like for the skew Poincaré observes.
Step14:  If Poincaré runs this experiment for a year, what is the chance he gets a statistically significant result?
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style='white')
from thinkstats2 import Pmf, Cdf
import thinkstats2
import thinkplot
decorate = thinkplot.config
Explanation: Poincaré's bread
Allen Downey
MIT License
End of explanation
def choose_bread(mean, std, n):
    loaves = np.random.normal(mean, std, size=n)
    return max(loaves)
mean = 950
std = 50
n = 5
choose_bread(mean, std, n)
Explanation: Poincaré's bread
Henri Poincaré was a French mathematician who taught at the Sorbonne around 1900.The following anecdote about him is probably fabricated, but it makes an interesting probability problem.
Supposedly, Poincaré suspected that his local bakery was selling loaves of bread that were lighter than the advertised weight of 1 kg, so every day for a year he bought a loaf of bread, brought it home and weighed it. At the end of the year, he plotted the distribution of his measurements and showed that it fit a normal distribution with mean 950 g and standard deviation 50 g. He brought this evidence to the bread police, who gave the baker a warning.
For the next year, Poincaré continued the practice of weighing his bread every day. At the end of the year, he found that the average weight was 1,000 g, just as it should be, but again he complained to the bread police, and this time they fined the baker.
Why? Because the shape of the distribution was asymmetric. Unlike the normal distribution, it was skewed to the right, which is consistent with the hypothesis that the baker was still making 950 g loaves, but deliberately giving Poincaré the heavier ones.
Write a program that simulates a baker who chooses n loaves from a distribution with mean 950 g and standard deviation 50 g, and gives the heaviest one to Poincaré. What value of n yields a distribution with mean 1000 g? What is the standard deviation?
Compare this distribution to a normal distribution with the same mean and the same standard deviation. Is the difference in the shape of the distribution big enough to convince the bread police?
First, I'll make a function that simulates choosing n loaves and giving Poincaré the heaviest.
End of explanation
sample = [choose_bread(mean, std, n) for i in range(10000)]
np.mean(sample)
Explanation: Here's a sample from the distribution of bread weight if the baker chooses 5 loaves.
End of explanation
for n in range(1, 10):
    sample = [choose_bread(mean, std, n) for i in range(10000)]
    print(n, np.mean(sample))
Explanation: And here's how the observed mean depends on n.
End of explanation
n = 4
sample = [choose_bread(mean, std, n) for i in range(10000)]
mu = np.mean(sample)
sigma = np.std(sample)
mu, sigma
Explanation: Looks like n=4 yields a measured mean close to 1000 g.
Here's a large sample of bread weights with n=4 so we can estimate the mean and standard deviation.
End of explanation
from scipy.stats import norm
xs = sample.copy()
xs.sort()
ys = norm.cdf(xs, mu, sigma)
plt.plot(xs, ys, color='gray', label='Gaussian model')
cdf = Cdf(sample, label='Sample with n=4')
thinkplot.Plot(cdf)
decorate(xlabel='Weight in g', ylabel='CDF', title='Distribution of bread weight')
Explanation: Now let's compare the distribution of the sample to the normal distribution with the same mean and std.
End of explanation
def simulate_year(mean, std, n, days):
    sample = [choose_bread(mean, std, n) for i in range(days)]
    return sample
Explanation: Visually they are almost indistinguishable, which makes it seem unlikely that Poincaré could really tell the difference, especially with only 365 observations.
But let's look more closely.  First, let's simulate one year.
End of explanation
from scipy.stats import skew
def test_stat(sample):
    return skew(sample)
Explanation: As a test statistic, we'll use the sample skew.
End of explanation
def run_model(mu, sigma, days):
    return np.random.normal(mu, sigma, size=days)
Explanation: The null hypothesis is that the data come from a normal distribution with the observed mean and standard deviation.
End of explanation
test_stat(run_model(mu, sigma, 365))
Explanation: Now we can generate fake data under the null hypothesis and compute the test statistic.
End of explanation
skew_under_null = [test_stat(run_model(mu, sigma, 365)) for i in range(1000)]
cdf = Cdf(skew_under_null)
thinkplot.Cdf(cdf)
decorate(xlabel='Skew', ylabel='CDF', title='Distribution of skew under H0')
Explanation: And we can estimate the distribution of the test statistic under the null hypothesis:
End of explanation
thresh = np.percentile(skew_under_null, 95)
thresh
Explanation: Here's the 95th percentile of the test statistic under the null hypothesis.
End of explanation
n = 4
skew_poincare = [test_stat(simulate_year(mean, std, n, 365)) for i in range(1000)]
Explanation: If Herni observes a skew above this threshold, he would conclude that it is statistically significant.
Let's see how often that would happen under the alternate hypothesis, that they baker is giving him the heaviest of 4 loaves.
End of explanation
cdf = Cdf(skew_poincare)
thinkplot.Cdf(cdf)
decorate(xlabel='Skew', ylabel='CDF', title='Distribution of skew if the baker cheats n=4')
Explanation: Here's what the sampling distribution looks like for the skew Poincaré observes.
End of explanation
np.mean(skew_poincare > thresh)
Explanation: If Poincaré runs this experiment for a year, what is the chance he gets a statistically significant result?
End of explanation
<END_TASK>