<jupyter_start><jupyter_text>Lambda School Data Science Module 133 Introduction to Bayesian Inference  Assignment - Code it up!Most of the above was pure math - now write Python code to reproduce the results! This is purposefully open ended - you'll have to think about how you should represent probabilities and events. You can and should look things up, and as a stretch goal - refactor your code into helpful reusable functions!Specific goals/targets:1. Write a function `def prob_drunk_given_positive(prob_drunk_prior, prob_positive, prob_positive_drunk)` that reproduces the example from lecture, and use it to calculate and visualize a range of situations2. Explore `scipy.stats.bayes_mvs` - read its documentation, and experiment with it on data you've tested in other ways earlier this week3. Create a visualization comparing the results of a Bayesian approach to a traditional/frequentist approach4. In your own words, summarize the difference between Bayesian and Frequentist statisticsIf you're unsure where to start, check out [this blog post of Bayes theorem with Python](https://dataconomy.com/2015/02/introduction-to-bayes-theorem-with-python/) - you could and should create something similar!Stretch goals:- Apply a Bayesian technique to a problem you previously worked (in an assignment or project work) on from a frequentist (standard) perspective- Check out [PyMC3](https://docs.pymc.io/) (note this goes beyond hypothesis tests into modeling) - read the guides and work through some examples- Take PyMC3 further - see if you can build something with it! ---  Resources - [Worked example of Bayes rule calculation](https://en.wikipedia.org/wiki/Bayes'_theoremExamples) (helpful as it fully breaks out the denominator)- [Source code for mvsdist in scipy](https://github.com/scipy/scipy/blob/90534919e139d2a81c24bf08341734ff41a3db12/scipy/stats/morestats.pyL139) ---<jupyter_code># Imports
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
from scipy import stats<jupyter_output><empty_output><jupyter_text>---  1. Write a function that reproduces the example from lecture, and use it to calculate and visualize a range of situations`def prob_drunk_given_positive(prob_drunk_prior, prob_positive, prob_positive_drunk)`<jupyter_code># First, go through one iteration by hand before writing into a function
prob_pos = 0.08
prob_pos_drunk = 1 - prob_pos
prob_drunk_prior = 1 / 1000
prob_NotThatDrunk = 1 - prob_drunk_prior<jupyter_output><empty_output><jupyter_text>$$P(Drunk|+) = \frac{prob\_pos\_drunk \cdot prob\_drunk\_prior}{(prob\_pos\_drunk \cdot prob\_drunk\_prior) + (prob\_positive \cdot prob\_NotThatDrunk)}$$<jupyter_code>(prob_pos_drunk * prob_drunk_prior) / ((prob_pos_drunk * prob_drunk_prior) + (prob_pos * (1 - prob_drunk_prior)))
# Now that I have everything laid out, writing the equation into a function
def prob_drunk_given_positive(prob_drunk_prior, prob_positive):
    return (((1 - prob_positive) * prob_drunk_prior) / 
            (((1 - prob_positive) * prob_drunk_prior) + (prob_positive * (1 - prob_drunk_prior))))
# Test out the function, one time!
prob_drunk_prior = prob_drunk_given_positive(prob_drunk_prior, prob_pos)
prob_drunk_prior
# Test out the function, ten times!
prior = 1 / 1000  # Reset the prior whenever cell is run

prior_data = []

for x in range(10):
    prior = prob_drunk_given_positive(prior, prob_pos)
    prior_data.append(prior)
    
print(prior_data)
# Graph the probabilities
sns.lineplot(data=pd.DataFrame(prior_data));
plt.title("Probability that driver is drunk given x positive breathalyzer tests");<jupyter_output><empty_output><jupyter_text>---  2. Explore `scipy.stats.bayes_mvs`- Read its [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.bayes_mvs.htmlscipy.stats.bayes_mvs)- Experiment with it on data you've tested in other ways earlier this week<jupyter_code>?stats.bayes_mvs<jupyter_output><empty_output><jupyter_text>ExperimentationalistFor this experimentation section, I'll use the `house-votes-84.data` dataset.<jupyter_code># Get the dataset from UCI
!wget https://archive.ics.uci.edu/ml/machine-learning-databases/voting-records/house-votes-84.data
# Load the dataset into a dataframe and add column headers
df = pd.read_csv('house-votes-84.data', 
                 header=None,
                 names=['party','handicapped-infants','water-project',
                          'budget','physician-fee-freeze', 'el-salvador-aid',
                          'religious-groups','anti-satellite-ban',
                          'aid-to-contras','mx-missile','immigration',
                          'synfuels', 'education', 'right-to-sue','crime','duty-free',
                          'south-africa'])

# Replace "?" with NaN and "y" / "n" with 1 / 0 
df = df.replace({
    "?": np.NaN,
    "n": 0,
    "y": 1,
})

print(df.shape)
df.head()<jupyter_output>(435, 17)<jupyter_text>> Calculate 95% confidence interval for "education" using Bayesian approach<jupyter_code># First, can't have any null values!
# Let's count those bad boiz
df["education"].isnull().value_counts()  
# 31 abstentions, resulting in 31 NaN values
# Some senators just don't care about the kids...
# Create pd.Series without the absentees
# Why yes...I do like to break things up into as many steps as possible (sometimes)
edumacation = df["education"].dropna(how="any")
print(len(edumacation))  # 404: government not found. Please try again later...
mean_ci_edu, _, _ = stats.bayes_mvs(edumacation, alpha=0.95)
print(mean_ci_edu)<jupyter_output>Mean(statistic=0.42326732673267325, minmax=(0.37488391760148726, 0.47165073586385925))<jupyter_text>> Calculate 99% confidence interval for "duty-free" also using Bayesian approach<jupyter_code># Now to do that all in one go, because I can
# ...and I like experimentating
mean_ci_doody, _, _ = stats.bayes_mvs(df["duty-free"].dropna(), alpha=0.99)
print(mean_ci_doody)<jupyter_output>Mean(statistic=0.4275184275184275, minmax=(0.3639768162425481, 0.49106003879430693))<jupyter_text>> Calculate 99% Bayesian confidence interval for "education", by party<jupyter_code># Create new dataframe with party as columns - so I can loop over to create the intervals
# Ready for some hot and steamy `.melt()` action?
# I don't even need to...doesn't really make sense in this instance.
print("I just wanted to say that...")
party_education = df[["party", "education"]].dropna()  # Also just wanted to say that
party_education.head()
print(type(df["duty-free"]))
# Iterate through the parties, calculating the Bayesian confidence intervals for each

confidence_level = 0.99

# Group by party, then iterate over the groups
for name, group in party_education.groupby(party_education["party"]):
    mean_ci_group, _, _ = stats.bayes_mvs(group["education"], alpha=confidence_level)
    print(f"Party: {name.title()}")
    print(f"Mean: {mean_ci_group[0]}")
    print(f"{confidence_level * 100}% Confidence Interval: {mean_ci_group[1]}")
    print("")<jupyter_output>Party: Democrat
Mean: 0.14457831325301204
99.0% Confidence Interval: (0.08661050929949306, 0.20254611720653104)

Party: Republican
Mean: 0.8709677419354839
99.0% Confidence Interval: (0.8005115457663251, 0.9414239381046426)<jupyter_text>---  3. Create a visualization comparing the results of a Bayesian approach to a traditional/frequentist approach<jupyter_code># Frequentist calculation of confidence interval

def confidence_interval(data, confidence=0.95):
    """
    Calculates a confidence interval around a sample mean for given data,
    using t-distribution and two-tailed test, default 95% confidence. 

    Arguments:
    data - iterable (list or numpy array) of sample observations
    confidence - level of confidence for the interval

    Returns:
    tuple of (mean, lower bound, upper bound)
    """
    data = np.array(data)  # Standardize to numpy array
    mean = np.mean(data)
    n = len(data)
    
    # Stdev divided by sqare root of n (degrees of freedom / # observations)
    stderr = np.std(data, ddof=1) / np.sqrt(n)
    # stderr = stats.sem(data)
    
    # Std error multiplied by t-statistic
    margin_of_error = stderr * stats.t.ppf((1 + confidence) / 2.0, n - 1)
    print(margin_of_error)

    return (mean, mean - margin_of_error, mean + margin_of_error)
# This time, I'm a Freq for confidence intervals
ci_freq = confidence_interval(edumacation)
print(ci_freq)
# I already hit the BAEsian, but tuple-ize it to plot it
ci_bae = tuple(mean_ci_edu)
print(ci_bae)
# Plot with seaborn kdeplot
sns.kdeplot(edumacation)

# Bayesian CI
plt.axvline(x=ci_bae[1][0], color="r");
plt.axvline(x=ci_bae[1][1], color="r");
plt.axvline(x=ci_bae[0], color="orange");

# Frequentist CI
plt.axvline(x=ci_freq[1], color="b");
plt.axvline(x=ci_freq[2], color="b");
plt.axvline(x=ci_freq[0], color="purple");
# Based on the plot given as an example in the `.bayes_mvs()` documentation

# Create figure and axes
fig = plt.figure()
ax = fig.add_subplot(111)

ax.hist(edumacation, bins=100, density=True, label='Histogram of data')

ax.vlines(ci_bae[0], 0, 0.5, colors='r', label='Estimated mean (Bayesian)')
ax.axvspan(ci_bae[1][0],ci_bae[1][1], facecolor='r',
           alpha=0.2, label=r'Estimated mean (99% limits)')

ax.vlines(ci_freq[0], 0, 0.5, colors='b', label='Estimated mean (Frequentist)')
ax.axvspan(ci_freq[1],ci_freq[2], facecolor='b',
           alpha=0.2, label=r'Estimated mean (99% limits)')

ax.legend(fontsize=10)
# ax.set_xlim([-4, 4])
# ax.set_ylim([0, 0.5])

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11241242287984941, 'fraction_numerical': 0.04005019345393705, 'mean_word_length': 4.778851963746224, 'pattern_counts': {'":': 3, '<': 33, '<?xml version=': 0, '>': 36, 'https://': 6, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '435713', 'license': 'MIT', 'n_tokens_mistral': 3487, 'n_tokens_neox': 3101, 'n_words': 904, 'path': 'module3-introduction-to-bayesian-inference/LS_DS_133_Introduction_to_Bayesian_Inference_Assignment.ipynb', 'repo_name': 'tobias-fyi/DS-Unit-1-Sprint-3-Statistical-Tests-and-Experiments', 'size': '9563'}