```python
from scipy.stats import norm
from scipy.stats import t
import numpy as np
import pandas as pd
from numpy.random import seed
import matplotlib.pyplot as plt
```


```python
#generates random variates to draw five samples from the standard normal distribution
sample1 = norm.rvs(size=5)
```


```python
np.mean(sample1)
```




    0.20789409052637584




```python
seed(47)
sample2 = norm.rvs(size=5)
sample2
```




    array([-0.84800948,  1.30590636,  0.92420797,  0.6404118 , -1.05473698])




```python
np.mean(sample2)
```




    0.19355593334131074




```python
#standard deviation - calculating manually
diff_squared = (sample2 - np.mean(sample2))**2
v = np.sum(diff_squared) / len(sample2)
print('Variance: ' + str(v))
print('Standard Deviation: ' + str(np.sqrt(v)))
```

    Variance: 0.9227899466393845
    Standard Deviation: 0.9606195639478641



```python
#standard deviation - calculating using np.sqrt(np.var())
print('Variance: ' + str(np.sqrt(np.var(sample2))))
```

    Variance: 0.9606195639478641



```python
#standard deviation - calculating using np.std()
print('Standard Deviation: ' + str(np.std(sample2)))
```

    Standard Deviation: 0.9606195639478641



```python
#standard deviation for sample
var_b = np.sum(diff_squared) / (len(sample2)-1)
np.sqrt(var_b)
```




    1.0740053227518152




```python
#with np.std() we get population standard deviation
np.std(sample2)
```




    0.9606195639478641




```python
#adding ddof=1 to np.std() will return the sample standard deviation
np.std(sample2, ddof=1)
```




    1.0740053227518152



For the sampling distribution of the mean, the standard deviation of this distribution is given by

\begin{equation}
\sigma_{mean} = \frac{\sigma}{\sqrt n}
\end{equation}

where $\sigma_{mean}$ is the standard deviation of the sampling distribution of the mean and $\sigma$ is the standard deviation of the population (the population parameter).

Let us imagine we live in a town of 50000 people and we know the height of everyone in this town. We will have 50000 numbers that tell us everything about our population. We'll simulate these numbers now and put ourselves in one particular town, called 'town 47', where the **population mean** height is **172** cm and **population standard deviation** is **5** cm.


```python
#Generate 50000 samples from the population heights using 172 population mean and 5 population std
seed(47)
pop_heights = norm.rvs(172, 5, size=50000)
pop_heights.min(), pop_heights.max(), np.mean(pop_heights), np.std(pop_heights)
```




    (150.86806415764312, 193.22214997313296, 172.0192602425845, 4.990839139566387)




```python
#Using histogram to show the distribution of heights in the population
_ = plt.hist(pop_heights, bins=30)
_ = plt.xlabel('height (cm)')
_ = plt.ylabel('number of people')
_ = plt.title('Distribution of heights in entire town population')
_ = plt.axvline(172, color='r')
_ = plt.axvline(172+5, color='r', linestyle='--')
_ = plt.axvline(172-5, color='r', linestyle='--')
_ = plt.axvline(172+10, color='r', linestyle='-.')
_ = plt.axvline(172-10, color='r', linestyle='-.')
```

Now, 50000 people is rather a lot to chase after with a tape measure. If all you want to know is the average height of the townsfolk, then can you just go out and measure a sample to get a pretty good estimate of the average height?


```python
#creating a function to create a sample size(n) from our population heights
def townsfolk_sampler(n):
    return np.random.choice(pop_heights, n)
```

Let's say you go out one day and randomly sample 10 people to measure.


```python
seed(47)
daily_sample1 = townsfolk_sampler(10)
```


```python
_ = plt.hist(daily_sample1, bins=10)
_ = plt.xlabel('height (cm)')
_ = plt.ylabel('number of people')
_ = plt.title('Distribution of heights in sample size 10')
```

The sample distribution doesn't look much like what we know (but wouldn't know in real-life) the population distribution looks like. What do we get for the mean?


```python
np.mean(daily_sample1)
```




    173.47911444163503



And if we went out and repeated this experiment?


```python
daily_sample2 = townsfolk_sampler(10)
```


```python
np.mean(daily_sample2)
```




    173.7317666636263



Simulate performing this random trial every day for a year, calculating the mean of each daily sample of 10, and plot the resultant sampling distribution of the mean.


```python
# Using a foor loop to perform this sampling of 10 for each day of a year
#[np.mean(townsfolk_sampler(10)) for x in range(365)]
trial_10 = []
for x in range(365):
    trial_10.append(np.mean(townsfolk_sampler(10)))
```


```python
seed(47)
# plotting the sampling distribution of the mean
_ = plt.hist(trial_10, bins = 10)
_ = plt.xlabel('height (cm)')
_ = plt.ylabel('number of people')
_ = plt.title('Distribution of heights of daily sample of 10 for a year')
_ = plt.axvline(172, color='r')
_ = plt.axvline(172+5, color='r', linestyle='--')
_ = plt.axvline(172-5, color='r', linestyle='--')
_ = plt.axvline(172+10, color='r', linestyle='-.')
_ = plt.axvline(172-10, color='r', linestyle='-.')

print(len(trial_10))
print(np.mean(trial_10))
print(np.median(trial_10))
```

The above is the distribution of the means of samples of size 10 taken from our population. The Central Limit Theorem tells us the expected mean of this distribution will be equal to the population mean, and standard deviation will be  𝜎/𝑛⎯⎯√ , which, in this case, should be approximately 1.58.

Q: Verify the above results from the CLT.


```python
np.std(trial_10, ddof=1)
```




    1.578160835870796



Remember, in this instance, we knew our population parameters, that the average height really is 172 cm and the standard deviation is 5 cm, and we see some of our daily estimates of the population mean were as low as around 168 and some as high as 176.

Q: Repeat the above year's worth of samples but for a sample size of 50 (perhaps you had a bigger budget for conducting surveys that year!) Would you expect your distribution of sample means to be wider (more variable) or narrower (more consistent)? Compare your resultant summary statistics to those predicted by the CLT.


```python
seed(47)
trial_50 = []
for x in range(365):
    trial_50.append(np.mean(townsfolk_sampler(50)))
```


```python
_ = plt.hist(trial_50, bins = 50, histtype='stepfilled', alpha=0.8)
_ = plt.xlabel('height (cm)')
_ = plt.ylabel('number of people')
_ = plt.title('Distribution of heights of daily sample of 50 for a year')
_ = plt.axvline(172, color='r')
_ = plt.axvline(172+5, color='r', linestyle='--')
_ = plt.axvline(172-5, color='r', linestyle='--')
_ = plt.axvline(172+10, color='r', linestyle='-.')
_ = plt.axvline(172-10, color='r', linestyle='-.')

print(len(trial_50))
print(np.mean(trial_50))
print(np.median(trial_50))
print(np.std(trial_50,ddof = 1))
```

What we've seen so far, then, is that we can estimate population parameters from a sample from the population, and that samples have their own distributions. Furthermore, the larger the sample size, the narrower are those sampling distributions.

### Let's now start from the position of knowing nothing about the heights of people in our town.
* Use our favorite random seed of 47, to randomly sample the heights of 50 townsfolk
* Estimate the population mean using np.mean
* Estimate the population standard deviation using np.std (remember which denominator to use!)
* Calculate the (95%) [margin of error](https://www.statisticshowto.datasciencecentral.com/probability-and-statistics/hypothesis-testing/margin-of-error/#WhatMofE) (use the exact critial z value to 2 decimal places - [look this up](https://www.statisticshowto.datasciencecentral.com/probability-and-statistics/find-critical-values/) or use norm.ppf())
* Calculate the 95% Confidence Interval of the mean
* Does this interval include the true population mean?


```python
seed(47)
second_trial_50 = townsfolk_sampler(50)
```


```python
np.mean(second_trial_50)
```




    172.7815108576788




```python
# Using ddof = 1 because this is a sample not population
np.std(second_trial_50, ddof=1)
```




    4.195424364433547




```python
critical_val = (norm.ppf(0.95))
# Standard Error of the Mean (a.k.a. the standard deviation of the sampling distribution of the sample mean!
se = (np.std(second_trial_50, ddof=1)) / (np.sqrt(50))
moe = critical_val * se # Margin of Error
print(moe)
```

    0.9759288364989565



```python
lower=np.mean(second_trial_50) - moe
upper=np.mean(second_trial_50) + moe
print('lower: ' + str(lower) + '\n' + 'Upper: ' + str(upper))
```

    lower: 171.80558202117984
    Upper: 173.75743969417775


__Q:__ Above we calculated the confidence interval using the critical z value. What is the problem with this? What requirement, or requirements, are we (strictly) failing?

__A:__ Because we don't know anything about our population, the dataset has an unknown population standard deviation. I think using a t-score instead of z-score would be better.

__Q__: Calculate the 95% confidence interval for the mean using the _t_ distribution. Is this wider or narrower than that based on the normal distribution above? If you're unsure, you may find this [resource](https://www.statisticshowto.datasciencecentral.com/probability-and-statistics/confidence-interval/) useful. For calculating the critical value, remember how you could calculate this for the normal distribution using norm.ppf().


```python
#Subtract 1 from your sample size to find the degrees of freedom (df).
dof = 50 - 1
t = t.ppf(0.95, (dof))
t
```




    1.6765508919142629




```python
# Margin of error for t-score
moe_t = t * se 
```


```python
lower=np.mean(second_trial_50) - moe_t
upper=np.mean(second_trial_50) + moe_t
print('lower: ' + str(lower) + '\n' + 'Upper: ' + str(upper))
```

    lower: 171.78677531740482
    Upper: 173.77624639795278


__A:__ Wider than the normal distribution since t-distribution curves are wider and bigger at the tails compared to normal distribution curves


```python

```

af2cf8bc2cee3e158885eef1139fb2192b81c0ec
ipynb
Jupyter Notebook
Python Statistics/Mini Projects/inferential_statistics_frequentist_mini-projects6.28.19/Mini Projects - Simplified.ipynb
atalebizadeh/DS-Career-Track
8bf78ef11041aef94810a392022cd51b94462d9c
MIT
Python Statistics/Mini Projects/inferential_statistics_frequentist_mini-projects6.28.19/Mini Projects - Simplified.ipynb
atalebizadeh/DS-Career-Track
8bf78ef11041aef94810a392022cd51b94462d9c
MIT
Python Statistics/Mini Projects/inferential_statistics_frequentist_mini-projects6.28.19/Mini Projects - Simplified.ipynb
atalebizadeh/DS-Career-Track
8bf78ef11041aef94810a392022cd51b94462d9c
MIT
Qwen/Qwen-72B
1. YES
2. YES