<jupyter_start><jupyter_text>Lab 7 - Inference<jupyter_code>% matplotlib inline<jupyter_output><empty_output><jupyter_text>This will make all the `matplotlib` images appear in the notebook.<jupyter_code>import numpy as np
import time
import seaborn as sns
import matplotlib.pyplot as plt
import scipy.stats as stats

sns.set(style="whitegrid")<jupyter_output><empty_output><jupyter_text>Directions**Failure to follow the directions will result in a "0"**The due dates for each are indicated in the Syllabus and the course calendar. If anything is unclear, please email <EMAIL> the official email for the course or ask questions in the Lab discussion area on Blackboard.The Labs also present technical material that augments the lectures and "book".  You should read through the entire lab at the start of each module. General Instructions1.  You will be submitting your assignment to Blackboard. If there are no accompanying files, you should submit *only* your notebook and it should be named using *only* your JHED id: fsmith79.ipynb for example if your JHED id were "fsmith79". If the assignment requires additional files, you should name the *folder/directory* your JHED id and put all items in that folder/directory, ZIP it up (only ZIP...no other compression), and submit it to Blackboard.        * do **not** use absolute paths in your notebooks. All resources should appear in the same directory as the rest of your assignments.    * the directory **must** be named your JHED id and **only** your JHED id.    2. Data Science is as much about what you write (communicating) as the code you execute (researching). In many places, you will be required to execute code and discuss both the purpose and the result. Additionally, Data Science is about reproducibility and transparency. This includes good communication with your team and possibly with yourself. Therefore, you must show **all** work.3. Avail yourself of the Markdown/Codecell nature of the notebook. If you don't know about Markdown, look it up. Your notebooks should not look like ransom notes. Don't make everything bold. Clearly indicate what question you are answering.4. Submit a cleanly executed notebook. The first code cell should say `In [1]` and each successive code cell should increase by 1 throughout the notebook.  Bayesian InferenceReally, there are only a few classical problems in statistical inference. However, we use the Bayes Theorem as the basis for solving all of them:$$P(H|D) = \frac{P(D|H)P(H)}{P(D)}$$You only need to identify what $H$ relates to...what is it? Is it some parameter of a distribution? Some property of a model (coefficients, error rate, etc.). For some formulations, we are more specific and specify $H$ as some parameter or parameters, $\theta$:$$P(\theta|D) = \frac{P(D|\theta)P(\theta)}{P(D)}$$In the text we saw how we could estimate the posterior distribution using four methods: Grid, Exact, Monte Carlo and Bootstrap. For this Lab, we'll concentrate on the Bootstrap method for the reasons specified in the text.  Statistical inference of a proportion in a Bernoulli Trial**1\. Suppose we have a coin that shows up heads 60% of the time ($\theta=p=0.6$). Generate 100 samples from this Binomial distribution (either as True/False or 1/0).**<jupyter_code>np.random.seed([1244875])

theta = 0.6
data = [1 if np.random.rand() < theta else 0 for _ in range( 100)]
print( data[0:20])<jupyter_output>[1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0]<jupyter_text>This is the synthetic data. At this point, we pretend that this is data we collected from the real world and we have no idea what $\theta$ really is. Understanding that inference is not certain, our goal is to make inferences about this parameter's value using this data we just "collected." Normally, the first thing we do is just calculate the parameter from our data. An *estimate* of some real world parameter is often given a "hat", for example $\theta$ becomes $\hat{\theta}$. Sometimes, it goes from Greek to Latin as in $\sigma$ to $s$ and sometimes it gets an adornment as well as in $\mu$ to $\bar{x}$.**2\. Calculate $\hat{theta}$.**<jupyter_code>theta_est = np.mean( data)
print( theta_est)<jupyter_output>0.67<jupyter_text>But we know that this $\hat{\theta}$ is not necessarily the "true" value. We want to know all the values that are supported by the data we collected and the degree to which they are supported...how confident we are in them. This is basically what we get when we calculate a posterior distribution over $\theta$ based on the data.And this is where the **(Non-Parameteric Bayesian) Bootstrap** estimate of that posterior distribution comes in. In the text we established *theoretically* how we went from a single data set to an estimate of the posterior distribution of our parameters. Now we're going to do it for reals.  Use the data we have to "bootstrap" an estimate of the posterior probability distribution over $\theta$, $P(\theta|D)$ which is "given the data we observed, how much are we to believe in the various values of $\theta$ and how much should we believe in them?" Remember that belief is quantified as probability.**3\. Generate the Bootstrap of the posterior distribution of $\hat{\theta}$ and answer the following questions:**First, we write a simple function to do our bootstrap sampling for us. It takes the data, a metric function and the number of bootstrap samples as the arguments. A metric function can be anything we like but it will most likely be something like `np.mean`, `np.var`, etc., it is whatever function we use to calculate our parameter/statistics.<jupyter_code>def bootstrap_sample( data, f, n=100):
    result = []
    m = len( data)
    for _ in range( n):
        sample = np.random.choice( data, len(data), replace=True)
        r = f( sample)
        result.append( r)
    return np.array( result)<jupyter_output><empty_output><jupyter_text>Now we used the function by supplying the data we "collected", our metric function `np.mean` and indicate we want 1000 bootstrap samples. This returns the data we can use as our posterior distribution of the proportion.<jupyter_code>posterior = bootstrap_sample( data, np.mean, 1000)<jupyter_output><empty_output><jupyter_text>If we like, we can plot a histogram of this posterior distribution:<jupyter_code>figure = plt.figure(figsize=(10, 6)) # first element is width, second is height.

axes = figure.add_subplot(1, 1, 1)

axes.hist( posterior, normed=True)
axes.set_ylabel( "Density")
axes.set_xlabel( "$\hat{theta}$")
axes.set_title( "Posterior Distribution of $\hat{theta}$")

plt.show()
plt.close()<jupyter_output><empty_output><jupyter_text>Note that while the data is discrete and boolean (true/false), the parameter $\theta$ is continuous. You might also notice that our distribution appears to be normally distributed. Based on the Central Limit Theorem, this is what we'd expect. **4\. What is the 90% Credible Interval (Bayesian Confidence Interval) for $\hat{\theta}$? Interpret it.**Although we'll often plot the posterior distribution, the real payoff from having it is to be able to do computations with it. There are a number of functions we can use for that purpose, for example, `mquantiles`. `mquantiles` is normally used to summarize the distributions of data but in this case, our data is estimates of $\theta$.<jupyter_code>stats.mstats.mquantiles( posterior, [0.05, 0.95])<jupyter_output><empty_output><jupyter_text>An important part of Data Science and assignments in this course is interpreting the results. This is not purely a coding class. Therefore, you should always, *always* interpret your results:There is a 90% probability that the value of $\theta$ is between 0.59 and 0.74 based on the data.Of course, there's nothing magical about only looking at the 90% confidence/credible interval and you can look at other ranges of interest as well. **5\. In Bayesian Statistics, we often identify a range of possible values for a parameter that we consider the same. This is known as the ROPE (Region of Practical Equivalance). We know that a fair coin would have $\theta$ of 0.5 but we're unlikely to get an exact value of 0.5 from our data. If the ROPE is 0.48-0.52, what is the probability that our coin's $\theta$ lies in that range and is thus "fair"?**<jupyter_code>np.mean((0.48 <= posterior) & (posterior <= 0.52))<jupyter_output><empty_output><jupyter_text>One of the downsides to the Boostrap approach is that we do not follow "Cromwell's Dictum" and we can get events with zero probability. We should just interpret these events are really have very small probabilities.Of course, now that we have this posterior distribution we can answer all kinds of (possibly) interesting and relevant questions to our problem. Let's stick with the basics, for now.  Exercises**Exercise 1.**In addition to estimates of the posterior distribution of parameters such as $\theta$, we are often interested in the posterior distribution of the *difference* of two $\theta$s. For example, we might be interested in the proportion of men who smoke versus the proportion of women who smoke.We can model $\theta_{men}$ and $\theta_{women}$ just as if they were coin biases and generate synthetic data just as if they were coin flips. Using the Non-Parametric Bootstrap, we can generate posterior distributions for $\hat{\theta}_{men}$ and $\hat{\theta}_{women}$ as well as $d$, the *difference*.These are the steps:1. Generate synthetic data for using $\theta_{men}$ = 0.23 and $\theta_{women}$ = 0.34,  with 100 observations each.2. Generate the bootstrap data for each.3. Generate difference data. You can do this by simply subtracting, element by element, one bootstrap sample from the other, $\theta_{men}$ - $\theta_{women}$.4. Plot the distributions of all three.5. Calculate the 90% Bayesian Confidence Interval of all three **and interpret them**.6. Determine a ROPE for the difference and tell me what's the probability that the "true" value of the difference falls in the ROPE.Use as many Markdown Cells and Code Cells as you need; it should look nice (not like a ransom note). 1. Generate random values between 0~12. For women:  If the value is <= 0.23, then 1 else 0 3. For men:  If the value is <= 0.34, then 1 else 0<jupyter_code>np.random.seed([111111])

theta = 0.23
women_data = [1 if np.random.rand() < theta else 0 for _ in range( 100)]

theta = 0.34
men_data = [1 if np.random.rand() < theta else 0 for _ in range( 100)]


print ("women: ",np.mean(women_data))
print ("men: ",np.mean(men_data))<jupyter_output>women:  0.23
men:  0.32<jupyter_text>The basic idea of using bootstrap sampling to estimate a posteriod distribution will stay with us throughout the entire semester. This will be our fundamental approach to statistical inference (there are other approaches and there are other *Bayesian* approaches). The important thing is to understand 1. why and 2. the dimensions along which the problems can vary such as,1. The nature of data. The data may take on a variety of different types. We've looked primarily at boolean or Bernoulli data. However, the data might be categorical (more than two discrete outcomes), counts, real valued, etc. This means that there may be more than one $\theta$. For example, the normal distribution has two $\theta$s: the mean, $\mu$, and the variance, $\sigma^2$. But you should think even more broadly than this. A linear regression as many $\theta$s: the coefficients, the coefficient of determination, the error of the regression, etc. A decision tree has a structure and error rate.2. A related concept is variability. We may have two true values, 0.23 and 0.24, but the variability of the data may not permit us to distinguish between them.3. Another dimension is the amount of data. We may not be able to get a "good" inference because we have not collected enough data.And, of course, all of these will and do interact. And a lot of experimental design is based on trying to limit variability (by "holding other things constant") and to get the "right" amount of data to support the inference we want to make.These exercises investigate some of the dimensions. **Exercise 2****1\. Repeat the guide example (coin flips) with a $\theta = 0.05$ and discuss. Were the credible intervals the same size? Was your estimate of $\theta$ as good? What does this say about statistical inference on relatively rare events or extreme values?**<jupyter_code>np.random.seed([87928356])<jupyter_output><empty_output><jupyter_text>** Statistical Inference for a single real valued $\theta$****Exercise 3**We can do the same thing for a real valued data (like weights, heights, etc.) and the $\theta$'s that describe such distributions. If we have a normal distribution, there are two such $\theta$s, $\mu$, the mean, and $\sigma$, the standard deviation. Remember, however, that we often think of the dispersion of our data as a percent of the mean or the *coefficient of variation*, v. **a\. Generate 50 observations from a normal distribution with $\mu=102.7$ and $v=5\%$.**You should refer to the previous Lab for generating synthetic data from the normal distribution and working with $v$, the coefficient of variation.<jupyter_code>np.random.seed([2386431651])<jupyter_output><empty_output><jupyter_text>**b. What is $\bar{x}$?** **c. Generate the Bootstrap estimate of the posterior distribution of $\bar{x}$.** **d. What is the 90% Credible Interval (Bayesian Confidence Interval) for $\hat{\theta}$?** **e. Define a ROPE of "about 100". What is the probability that $\bar{x}$ falls within the ROPE?.** **Exercise 4\. Repeat Steps 1-5 with $v=25\%$.**<jupyter_code>np.random.seed([484716248])<jupyter_output><empty_output><jupyter_text>**Exercise 5\. Repeat Steps 1-5 with $v=25\%$ and 500 samples.**<jupyter_code>np.random.seed([484716248])<jupyter_output><empty_output><jupyter_text>** Statistical Inference for a two real valued $\theta$s****Exercise 6. Following the text, apply the Bootstrap to make inferences about the difference between $\mu_1$ and $\mu_2$**1. Data set 1 has $\mu_1=102.7$ and $v_1=10\%$ and 100 observations.2. Data set 2 has $\mu_2=104.2$ and $v_2=5\%$ and 100 observations.Pay special consideration to formulating your ROPE for the *difference* between the two parameters and making inferences about it.<jupyter_code>np.random.seed([67366372])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08482825754415241, 'fraction_numerical': 0.02037268808232513, 'mean_word_length': 4.771669341894061, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 59, 'https://': 0, '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': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '69337', 'license': 'MIT', 'n_tokens_mistral': 4110, 'n_tokens_neox': 3768, 'n_words': 2055, 'path': 'lab7/Lab 7.ipynb', 'repo_name': 'topseer/APL_Great', 'size': '14394'}