<jupyter_start><jupyter_text>Download this page as a jupyter notebook at [Lesson 18](http://54.243.252.9/engr-1330-webroot/1-Lessons/Lesson18/ENGR-1330-Lesson18.ipynb)---<jupyter_code># Preamble script block to identify host, user, and kernel
import sys
! hostname
! whoami
print(sys.executable)
print(sys.version)
print(sys.version_info)
%%html
<!--Script block to left align Markdown Tables-->
<style>
  table {margin-left: 0 !important;}
</style><jupyter_output><empty_output><jupyter_text>--- ENGR 1330 Computational Thinking with Data Science Last GitHub Commit Date: 1 Nov 2021 Lesson 18 : Interval Estimates (Confidence and Prediction Intervals)A procedure to quantify the probability of an estimate being representative of an equivalent population value Objectives- To apply fundamental concepts involved in probability estimation modeling and descriptive statistics;  - Is a descriptive statistic a good estimate of its population value?  - What is the best estimate of a statistic from a set, and its likely range?--- Computational Thinking ConceptsThe CT concepts include:- Abstraction => Represent data behavior with a statistic - Pattern Recognition => Compare patterns in (our) data models to make a decision--- Textbook Resources[https://www.inferentialthinking.com/chapters/11/Testing_Hypotheses.html](https://www.inferentialthinking.com/chapters/11/Testing_Hypotheses.html) Much of what follows below is borrowed with some changes for our class from: https://pubs.usgs.gov/tm/04/a03/tm4a3.pdf Interval EstimatesSeveral summary statistics were presented that described key attributes of a dataset, including sample estimates such as $\bar x$ and $s^2$, of true and unknown population parameters, such as $\mu$, the population mean, and $\sigma^2$, the population variance.What follows are methods to describe the uncertainty or reliability of sample estimates. As an alternative to reporting a single estimate, the utility of reporting a range of values called an interval estimate is demonstrated. Both parametric and nonparametric interval estimates are presented. These intervals can also be used to test whether the population parameter is significantly different from some prespecified value. Definition of Interval EstimatesThe sample median and sample mean are two types of estimates of the central tendency of a population. Such estimates are called point estimates. By themselves, point estimates do not portray the reliability, or lack of reliability (variability), of these estimates. For example, suppose that two datasets—X and Y—exist, both have a sample mean of about 50 and contain the same number of observations.<jupyter_code>import numpy
import pandas
process1 = []
process2 = []
for i in range(1,100,1):
    process1.append(i)
for i in range(41,61,1):
    process2.append(i)
    
x = numpy.random.choice(process1,500)
y = numpy.random.choice(process2,500)

# make lists into dataframe 
boxdf = pandas.DataFrame()
boxdf['X']=x
boxdf['Y']=y

# describe lists/dataframe 
print(boxdf.describe())

boxdf.plot.hist(alpha=0.5); #<jupyter_output>X           Y
count  500.00000  500.000000
mean    49.81800   49.888000
std     26.93756    5.846008
min      1.00000   41.000000
25%     26.75000   45.000000
50%     50.00000   49.000000
75%     71.00000   55.000000
max     99.00000   60.000000<jupyter_text>The Y data all cluster tightly around 50; the X data are much more variable. The point estimate of 50 for dataset X is less certain than the point estimate for dataset Y because of the greater variability in the X data. Reporting only the sample (point) estimate of 50 fails to give any hint of this difference, although it is graphically quite evident.We would then recall we have been hypothesis testing the crap out of things, so we can apply those tools:<jupyter_code># are the samples approximatley normal?
# the Shapiro-Wilk Normality Test for X
from scipy.stats import shapiro
stat, p = shapiro(boxdf['X'])
print('stat=%.3f, p=%.3f' % (stat, p))
if p > 0.05:
        print('Probably Normal')
else:
        print('Probably Abby Normal')
# are the samples approximatley normal?
# the Shapiro-Wilk Normality Test for Y
from scipy.stats import shapiro
stat, p = shapiro(boxdf['Y'])
print('stat=%.3f, p=%.3f' % (stat, p))
if p > 0.05:
        print('Probably Normal')
else:
        print('Probably Abby Normal')
from scipy.stats import mannwhitneyu # import a useful non-parametric test
stat, p = mannwhitneyu(boxdf['X'],boxdf['Y'])
print('statistic=%.3f, p-value at rejection =%.3f' % (stat, p))
if p > 0.05:
    print('Probably the same mean')
else:
    print('Probably different mean')<jupyter_output>statistic=124837.000, p-value at rejection =0.486
Probably the same mean<jupyter_text>At this point we know that the mean values are about the same, but the two datasets are not alike (its the variance stupid!)As an alternative to just point estimates, interval estimates are intervals that have a stated probability of containing the true population value. In our example, the mean for Y is 50 and we are pretty sure of that, because the sample variance is small - if our sample had also tested as normal we could state that the estimate of population mean is 50 +/- 18 ( 3 standard deviations) with probability of 99%. Such a statement is an interval estimate.In general, we will be presenting two-sided intervals (where the probability of the true value being higher than the upper limit is assumed to be equal to the probability of its being lower than the lower limit). There are also one-sided intervals. The intervals are wider for datasets having greater variability and the same number of data points.Thus, in the above example, an interval between 38 and 62 may have a 95-percent probability of containing the (unknown) true population mean of dataset Y. It would take a much wider interval, say between -6 and 106, to have the same probability of containing the true mean of dataset X.The difference in the reliability of the two estimates is therefore clearly stated using interval estimates. Interval estimates can provide two pieces of information which point estimates cannot:   > A statement of the probability or likelihood that the interval contains the true population value (its reliability), and    > A statement of the likelihood that a single data point with specified magnitude comes from the population under study.Interval estimates for the first purpose are called confidence intervals; intervals for the second purpose are called prediction intervals. Though related, the two types of interval estimates are not identical, and should be interchanged. ---**First a visit with a textbook!**Much of what follows below is borrowed with minor changes for our class from: https://www.inferentialthinking.com/chapters/13/1/Percentiles.html PercentilesNumerical data can be sorted in increasing or decreasing order. Thus the values of a numerical data set have a rank order. A percentile is the value at a particular rank. We have exploited this is our plotting position formulas earlier.For example, if your score on a test is on the 95th percentile, a common interpretation is that only 5% of the scores were higher than yours. The median is the 50th percentile; it is commonly assumed that 50% the values in a data set are above the median.But some care is required in giving percentiles a precise definition that works for all ranks and all lists. To see why, consider an extreme example where all the students in a class score 75 on a test. Then 75 is a natural candidate for the median, but it's not true that 50% of the scores are above 75. Also, 75 is an equally natural candidate for the 95th percentile or the 25th or any other percentile. Ties – that is, equal data values – have to be taken into account when defining percentiles.You also have to be careful about exactly how far up the list to go when the relevant index isn't clear. For example, what should be the 87th percentile of a collection of 10 values? The 8th value of the sorted collection, or the 9th, or somewhere in between? Again we have already addressed this issue with our plotting position formulas earlier.A Numerical ExampleBefore giving a general definition of all percentiles, we will define the 80th percentile of a collection of values to be the smallest value in the collection that is at least as large as 80% of all of the values.For example, let's consider the sizes of the five largest continents – Africa, Antarctica, Asia, North America, and South America – rounded to the nearest million square miles.<jupyter_code>import datascience
sizes = datascience.make_array(12, 17, 6, 9, 7)<jupyter_output><empty_output><jupyter_text>The 80th percentile is the smallest value that is at least as large as 80% of the elements of sizes, that is, four-fifths of the five elements. That's 12:<jupyter_code>numpy.sort(sizes)<jupyter_output><empty_output><jupyter_text>The 80th percentile is a value on the list, namely 12. You can see that 80% of the values are less than or equal to it, and that it is the smallest value on the list for which this is true.Analogously, the 70th percentile is the smallest value in the collection that is at least as large as 70% of the elements of sizes. Now 70% of 5 elements is "3.5 elements", so the 70th percentile is the 4th element on the list. That's 12, the same as the 80th percentile for these data. Not terribly useful yet but wait there's more!The percentile function The `datascience.percentile` functionThe percentile function takes two arguments: a rank between 0 and 100, and a array. It returns the corresponding percentile of the array.<jupyter_code>datascience.percentile(70, sizes)<jupyter_output><empty_output><jupyter_text>--- The General Definition **Below is from the Data8 Textbook -- compare to our future lab on plotting position formulas**Let 𝑝 be a number between 0 and 100. The 𝑝-th percentile of a collection is the smallest value in the collection that is at least as large as p% of all the values.By this definition, any percentile between 0 and 100 can be computed for any collection of values, and it is always an element of the collection.In practical terms, suppose there are 𝑛 elements in the collection. To find the 𝑝-th percentile:- Sort the collection in increasing order.- Find p% of n: (𝑝/100)×𝑛 Call that 𝑘- If 𝑘 is an integer, take the 𝑘-th element of the sorted collection.- If 𝑘 is not an integer, round it up to the next integer, and take that element of the sorted collection. Example The table scores_and_sections contains one row for each student in a class of 359 students. The columns are the student's discussion section and midterm score.Path to file is `http://54.243.252.9/engr-1330-webroot/4-Databases/`<jupyter_code>import requests
remote_url="http://54.243.252.9/engr-1330-webroot/4-Databases/scores_by_section.csv"  # set the url
rget = requests.get(remote_url, allow_redirects=True)  # get the remote resource, follow imbedded links
open('scores_by_section.csv','wb').write(rget.content); # extract from the remote the contents, assign to a local file same name<jupyter_output><empty_output><jupyter_text>The table scores_and_sections contains one row for each student in a class of 359 students. The columns are the student's discussion section and midterm score.<jupyter_code>scores_and_sections = datascience.Table.read_table('scores_by_section.csv')
scores_and_sections<jupyter_output><empty_output><jupyter_text>If we make a histogram (using the Data8 tools from our textbook) can we approximate the 85%?  (just by looking it seems to be at 22-ish from the plot)<jupyter_code>scores_and_sections.select('Midterm').hist(bins=numpy.arange(-0.5, 25.6, 1))
datascience.percentile(85, scores_and_sections.column(1))<jupyter_output><empty_output><jupyter_text>So the eyeball value was not too bad! According to the percentile function, the 85th percentile was 22. To check that this is consistent with the definition, let's apply the definition directly.First, put the scores in increasing order:<jupyter_code>sorted_scores = numpy.sort(scores_and_sections.column(1))
len(sorted_scores)<jupyter_output><empty_output><jupyter_text>There are 359 scores in the array. So next, find 85% of 359, which is 305.15.<jupyter_code>0.85 * 359<jupyter_output><empty_output><jupyter_text>That's not an integer. By the definition, the 85th percentile is the 306th element of sorted_scores, which, by Python's indexing convention, is item 305 of the array.<jupyter_code># The 306th element of the sorted array
sorted_scores.item(305)<jupyter_output><empty_output><jupyter_text>That's the same as the answer we got by using percentile. In future, we will just use percentile. Quartiles The first quartile of a numercial collection is the 25th percentile. The terminology arises from the first quarter. The second quartile is the median, and the third quartile is the 75th percentile.For our scores data, those values are:<jupyter_code>print(datascience.percentile(25, sorted_scores))
print(datascience.percentile(50, sorted_scores))
print(datascience.percentile(75, sorted_scores))
pandas.DataFrame(sorted_scores).describe()<jupyter_output><empty_output><jupyter_text>--- Bootstrapping (Hugely Important, but has a steep jargon curve)A data scientist is using the data in a random sample to estimate an unknown parameter. She uses the sample to calculate the value of a statistic that she will use as her estimate.Once she has calculated the observed value of her statistic, she could just present it as her estimate and go on her merry way. But she knows that her random sample is just one of numerous possible random samples, and thus her estimate is just one of numerous plausible estimates.By how much could those estimates vary? To answer this, it appears as though she needs to draw another sample from the population, and compute a new estimate based on the new sample. But she doesn't have the cash to go back to the population and draw another sample.It looks as though the data scientist is stuck, but a method called the bootstrap (or just bootstrapping) can help out. Because it is not feasible to generate new samples from the population, the bootstrap generates new random samples by a method called resampling: the new samples are drawn at random from the original sample.  (Essentially, the orginal sample is now treated as a complete population - its a hack, but surprisingly useful)  Employee Compensation in the City of San FranciscoSF OpenData is a website where the City and County of San Francisco make some of their data publicly available. One of the data sets contains compensation data for employees of the City. These include medical professionals at City-run hospitals, police officers, fire fighters, transportation workers, elected officials, and all other employees of the City.<jupyter_code>remote_url="http://54.243.252.9/engr-1330-webroot/4-Databases/san_francisco_2015.csv"  # set the url
rget = requests.get(remote_url, allow_redirects=True)  # get the remote resource, follow imbedded links
open('san_francisco_2015.csv','wb').write(rget.content); # extract from the remote the contents, assign<jupyter_output><empty_output><jupyter_text>Compensation data for the calendar year 2015 are in the table sf2015.<jupyter_code>sf2015 = datascience.Table.read_table('san_francisco_2015.csv')
sf2015<jupyter_output><empty_output><jupyter_text>There is one row for each of 42,979 employees. There are numerous columns containing information about City departmental affiliation and details of the different parts of the employee's compensation package. Here is the row correspoding to the late Edward Lee, the Mayor at that time.<jupyter_code>sf2015.where('Job', datascience.are.equal_to('Mayor'))<jupyter_output><empty_output><jupyter_text>We are going to study the final column, Total Compensation. That's the employee's salary plus the City's contribution towards his/her retirement and benefit plans.Financial packages in a calendar year can sometimes be hard to understand as they depend on the date of hire, whether the employee is changing jobs within the City, and so on. For example, the lowest values in the Total Compensation column look a little strange.<jupyter_code>sf2015.sort('Total Compensation')<jupyter_output><empty_output><jupyter_text>For clarity of comparison, we will focus our attention on those who had at least the equivalent of a half-time job for the whole year. At a minimum wage of about \\$10 per hour, and 20 hours per week for 52 weeks, that's a salary of about \\$10,000.<jupyter_code>sf2015 = sf2015.where('Salaries', datascience.are.above(10000))

sf2015.num_rows<jupyter_output><empty_output><jupyter_text>Population and ParameterLet this table of just over 36,500 rows be our population. Here is a histogram of the total compensations.<jupyter_code>sf_bins = numpy.arange(0, 700000, 25000)
sf2015.select('Total Compensation').hist(bins=sf_bins)<jupyter_output>/opt/jupyterhub/lib/python3.8/site-packages/datascience/tables.py:5206: UserWarning: FixedFormatter should only be used together with FixedLocator
  axis.set_xticklabels(ticks, rotation='vertical')<jupyter_text>While most of the values are below \\$300,000, a few are quite a bit higher. For example, the total compensation of the Chief Investment Officer was almost \\$650,000. That is why the horizontal axis stretches to \\$700,000.<jupyter_code>sf2015.sort('Total Compensation', descending=True).show(2)<jupyter_output><empty_output><jupyter_text>Now let the parameter be the median of the total compensations.Since we have the luxury of having all of the data from the population, we can simply calculate the parameter:<jupyter_code>pop_median = datascience.percentile(50, sf2015.column('Total Compensation'))
pop_median<jupyter_output><empty_output><jupyter_text>The median total compensation of all employees was just over \\$110,300.From a practical perspective, there is no reason for us to draw a sample to estimate this parameter since we simply know its value. But in this section we are going to pretend we don't know the value, and see how well we can estimate it based on a random sample. A Random Sample and an EstimateLet us draw a sample of 500 employees at random without replacement, and let the median total compensation of the sampled employees serve as our estimate of the parameter.<jupyter_code>our_sample = sf2015.sample(500, with_replacement=False)
our_sample.select('Total Compensation').hist(bins=sf_bins)
est_median = datascience.percentile(50, our_sample.column('Total Compensation'))
print('Estimated Median = ',est_median)<jupyter_output>Estimated Median =  110566.68<jupyter_text>The sample size is large, so that the distribution of the sample approximates that of the population, and consequently the sample median is not very far from the population median (though of course it is not exactly the same).So now we have one estimate of the parameter. But had the sample come out differently, the estimate would have had a different value. We would like to be able to quantify the amount by which the estimate could vary across samples. That measure of variability will help us measure how accurately we can estimate the parameter.To see how different the estimate would be if the sample had come out differently, we could just draw another sample from the population, but that would be cheating. We are trying to mimic real life, in which we won't have all the population data at hand.Somehow, we have to get another random sample without sampling from the population. The Bootstrap: Resampling from the SampleWhat we do have is a large random sample from the population. As we know, a large random sample is likely to resemble the population from which it is drawn. This observation allows data scientists to lift themselves up by their own bootstraps: the sampling procedure can be replicated by sampling from the sample.Here are the steps of the bootstrap method for generating another random sample that resembles the population:  > - Treat the original sample as if it were the population.  > - Draw from the sample, at random with replacement, the same number of times as the original sample size.It is important to resample the same number of times as the original sample size. The reason is that the variability of an estimate depends on the size of the sample. Since our original sample consisted of 500 employees, our sample median was based on 500 values. To see how different the sample could have been, we have to compare it to the median of other samples of size 500.If we drew 500 times at random without replacement from our sample of size 500, we would just get the same sample back. By drawing with replacement, we create the possibility for the new samples to be different from the original, because some employees might be drawn more than once and others not at all.Why is this a good idea? Because the distribution of the original sample is likely to resemble the population (if we took enough instances), and the distributions of all the "resamples" are likely to resemble the original sample. Then by an appeal to the transitive axiom, the distributions of all the resamples are likely to resemble the population as well. A Resampled MedianRecall that when the sample method is used without specifying a sample size, by default the sample size equals the number of rows of the table from which the sample is drawn. That's perfect for the bootstrap! Here is one new sample drawn from the original sample, and the corresponding sample median.<jupyter_code>resample_1 = our_sample.sample()
resample_1.select('Total Compensation').hist(bins=sf_bins)
resampled_median_1 = datascience.percentile(50, resample_1.column('Total Compensation'))
resampled_median_1<jupyter_output><empty_output><jupyter_text>By resampling, we have another estimate of the population median. By resampling again and again, we will get many such estimates, and hence an empirical distribution of the estimates.<jupyter_code>resample_2 = our_sample.sample()
resampled_median_2 = datascience.percentile(50, resample_2.column('Total Compensation'))
resampled_median_2<jupyter_output><empty_output><jupyter_text>Bootstrap Empirical Distribution of the Sample MedianLet us define a function `bootstrap_median` that takes our original sample, the label of the column containing the variable, and the number of bootstrap samples we want to take, and returns an array of the corresponding resampled medians.Each time we resample and find the median, we replicate the bootstrap process. So the number of bootstrap samples will be called the number of replications.<jupyter_code>def bootstrap_median(original_sample, label, replications):
    """Returns an array of bootstrapped sample medians:
    original_sample: table containing the original sample
    label: label of column containing the variable
    replications: number of bootstrap samples
    modified by t.g. cleveland for missing dependencies from 
    original UCB script
    """
    just_one_column = original_sample.select(label)
    medians = datascience.make_array()
    for i in numpy.arange(replications):
        bootstrap_sample = just_one_column.sample()
        resampled_median = datascience.percentile(50, bootstrap_sample.column(0))
        medians = numpy.append(medians, resampled_median)
        
    return medians<jupyter_output><empty_output><jupyter_text>We now replicate the bootstrap process 5,000 times. The array bstrap_medians contains the medians of all 5,000 bootstrap samples. Notice that the code takes longer to run than our previous code. It has a lot of resampling to do!<jupyter_code>bstrap_medians = bootstrap_median(our_sample, 'Total Compensation', 5000) # change the 5 to 5000 for production run!<jupyter_output><empty_output><jupyter_text>Here is the histogram of the 5000 medians. The red dot is the population parameter: it is the median of the entire population, which we happen to know but did not use in the bootstrap process.<jupyter_code>resampled_medians = datascience.Table().with_column('Bootstrap Sample Median', bstrap_medians)

#median_bins=np.arange(100000, 130000, 2500)
#resampled_medians.hist(bins = median_bins)
resampled_medians.hist()

import matplotlib.pyplot 
matplotlib.pyplot.scatter(pop_median, 0,color='red', s=120 );<jupyter_output><empty_output><jupyter_text>It is important to remember that the red dot is fixed: it is \\$110,305.79, the population median. The empirical histogram is the result of random draws, and will be situated randomly relative to the red dot.Remember also that the point of all these computations is to estimate the population median, which is the red dot. Our estimates are all the randomly generated sampled medians whose histogram you see above. We want those estimates to contain the parameter – it they don't, then they are off.--- Do the Estimates Capture the Parameter?How often does the empirical histogram of the resampled medians sit firmly over the red dot, and not just brush the dot with its tails? To answer this, we must define "sit firmly". Let's take that to mean "the middle 95% of the resampled medians contains the red dot".Here are the two ends of the "middle 95%" interval of resampled medians:<jupyter_code>left = datascience.percentile(2.5, bstrap_medians)
right = datascience.percentile(97.5, bstrap_medians)
print("lower/upper bounds of 95-percent of sample range \n",left," ",right)<jupyter_output>lower/upper bounds of 95-percent of sample range 
 104877.51   114055.55<jupyter_text>The population median of \\$110,305 is indeed between these two numbers. The interval and the population median are shown on the histogram below.<jupyter_code>#median_bins=np.arange(100000, 130000, 2500)
#resampled_medians.hist(bins = median_bins)
resampled_medians.hist()

matplotlib.pyplot.plot(datascience.make_array(left, right), datascience.make_array(0, 0), color='yellow', lw=9, zorder=1)
matplotlib.pyplot.scatter(pop_median, 0, color='red', s=200, zorder=2);<jupyter_output>/opt/jupyterhub/lib/python3.8/site-packages/datascience/tables.py:5206: UserWarning: FixedFormatter should only be used together with FixedLocator
  axis.set_xticklabels(ticks, rotation='vertical')<jupyter_text>The "middle 95%" interval of estimates captured the parameter in our example. But was that a fluke?To see how frequently the interval contains the parameter, we have to run the entire process over and over again. Specifically, we will repeat the following process 100 times:> - Draw an original sample of size 500 from the population.> - Carry out 5,000 replications of the bootstrap process and generate the "middle 95%" interval of resampled medians.We will end up with 100 intervals, and count how many of them contain the population median.Spoiler alert: The statistical theory of the bootstrap says that the number should be around 95. It may be in the low 90s or high 90s, but not much farther off 95 than that.<jupyter_code># THE BIG SIMULATION: This one takes several minutes.
# Generate 100 intervals, in the table intervals
# modified by tgc for UCB dependencies not loaded as part of textbook
# test 26oct2020 on AWS server runtime == about 4 minutes
left_ends = datascience.make_array()
right_ends = datascience.make_array()

total_comps = sf2015.select('Total Compensation')

for i in numpy.arange(100):
    first_sample = total_comps.sample(500, with_replacement=False)
    medians = bootstrap_median(first_sample, 'Total Compensation', 5000)
    left_ends = numpy.append(left_ends, datascience.percentile(2.5, medians))
    right_ends = numpy.append(right_ends, datascience.percentile(97.5, medians))

intervals = datascience.Table().with_columns(
    'Left', left_ends,
    'Right', right_ends
)<jupyter_output><empty_output><jupyter_text>For each of the 100 replications, we get one interval of estimates of the median.<jupyter_code>intervals<jupyter_output><empty_output><jupyter_text>The good intervals are those that contain the parameter we are trying to estimate. Typically the parameter is unknown, but in this section we happen to know what the parameter is.<jupyter_code>pop_median<jupyter_output><empty_output><jupyter_text>How many of the 100 intervals contain the population median? That's the number of intervals where the left end is below the population median and the right end is above.<jupyter_code>intervals.where('Left', datascience.are.below(pop_median)).where('Right', datascience.are.above(pop_median)).num_rows<jupyter_output><empty_output><jupyter_text>It takes a few minutes to construct all the intervals, but try it again if you have the patience. Most likely, about 95 of the 100 intervals will be good ones: they will contain the parameter.It's hard to show you all the intervals on the horizontal axis as they have large overlaps – after all, they are all trying to estimate the same parameter. The graphic below shows each interval on the same axes by stacking them vertically. The vertical axis is simply the number of the replication from which the interval was generated.The red line is where the parameter is. Good intervals cover the parameter; there are about 95 of these, typically.If an interval doesn't cover the parameter, it's a dud. The duds are the ones where you can see "daylight" around the red line. There are very few of them – about 5, typically – but they do happen.Any method based on sampling has the possibility of being off. The beauty of methods based on random sampling is that we can quantify how often they are likely to be off.<jupyter_code>replication_number = numpy.ndarray.astype(numpy.arange(1, 101), str)
intervals2 = datascience.Table(replication_number).with_rows(datascience.make_array(left_ends, right_ends))

matplotlib.pyplot.figure(figsize=(8,8))
for i in numpy.arange(100):
    ends = intervals2.column(i)
    matplotlib.pyplot.plot(ends, datascience.make_array(i+1, i+1), color='gold')
matplotlib.pyplot.plot(datascience.make_array(pop_median, pop_median), datascience.make_array(0, 100), color='red', lw=2)
matplotlib.pyplot.xlabel('Median (dollars)')
matplotlib.pyplot.ylabel('Replication')
matplotlib.pyplot.title('Population Median and Intervals of Estimates');<jupyter_output><empty_output><jupyter_text>To summarize what the simulation shows, suppose you are estimating the population median by the following process:> - Draw a large random sample from the population.> - Bootstrap your random sample and get an estimate from the new random sample.> - Repeat the above step thousands of times, and get thousands of estimates.> - Pick off the "middle 95%" interval of all the estimates.That gives you one interval of estimates. Now if you repeat the entire process 100 times, ending up with 100 intervals, then about 95 of those 100 intervals will contain the population parameter.In other words, this process of estimation captures the parameter about 95% of the time.You can replace 95% by a different value, as long as it's not 100. Suppose you replace 95% by 80% and keep the sample size fixed at 500. Then your intervals of estimates will be shorter than those we simulated here, because the "middle 80%" is a smaller range than the "middle 95%". Only about 80% of your intervals will contain the parameter.Bootstrapping is often used to construct statistical tests (from known distributions). ------ Confidence Intervals for the MedianA confidence interval for the true population median may be computed in two ways: (1) without assuming the data follow any specific distribution (nonparametric), or (2) assuming they follow a distribution such as the lognormal (parametric). --- Nonparametric Interval Estimate for the MedianWe will consider two nonparametric approaches to interval estimates for the median. The first is based on the binomial distribution and the second is a bootstrap method, which is a general approach to many estimation problems. Neither approach requires assumptions about the distribution of the random variable.For the binomial approach we start by selecting the desired significance level α, which is the acceptable risk of not including the true median. One-half of this risk (α/2) is assigned to each end of the interval. To compute the confidence interval for the median we use the cumulative distribution function (cdf) of the binomial distribution. Learn about the cdf at https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binom.htmlHere's a way to learn, lets code up a script to plot the CDF, and see if it looks like a CDF, in this case we will plot the accumulated probability on the x-axis, and the random variable value on the y-axis.<jupyter_code>import numpy
import pandas
import scipy.stats 
import matplotlib.pyplot

n = 25
p = 0.5
prob = []
xval = []
for rv in range(0,n):
    x=float(rv)
    prob.append(scipy.stats.binom.cdf(x, n, p))
    xval.append(x)

matplotlib.pyplot.figure(figsize=(8,8))
matplotlib.pyplot.plot(xval, prob,  lw=5)
plottitle = 'Binomial CDF for n= ' + str(round(n,3))+ ' p = ' + str(round(p,3))
matplotlib.pyplot.title(plottitle)
matplotlib.pyplot.xlabel('Accumulated Probability')
matplotlib.pyplot.ylabel('Random Variable Value')
matplotlib.pyplot.show()<jupyter_output><empty_output><jupyter_text>To determine the 100 · (1 − α)-percent confidence interval we use the binomial distribution function to determine the critical values of the ranks of the dataset that correspond to α/2 and 1 − (α/2) points on the binomial cumulative distribution function. These critical values of the ranks are then associated with their data values to form the upper and lower confidence limits for the median. We then use the binomial distribution to answer the following question: How likely is it that the true population median, $c_{0.5}$, would be such that k of the n observed data would be above $c_{0.5}$ and n−k below $c_{0.5}$, where for example, k could be 0, 1, 2, 3,…, 25 out of n = 25? The binomial distribution with prob = 0.5 is used because the probability of being above the population median is 0.5. The resulting confidence interval will reflect the shape (skewed or symmetric) of the original data. Nonparametric intervals cannot exactly produce the desired confidence level when the sample sizes are small. This is because the possible values are discrete, jumping from one data value to the next at the ends of the intervals. However, confidence intervals close to those desired are available for all but the smallest sample sizes. The process of computing the confidence interval is best illustrated with an example. Example: 25 arsenic concentrations were reported for groundwaters of southeastern New Hampshire (Boudette and others, 1985). A boxplot of the data is shown. Compute the α = 0.05 interval estimate of the median concentration.<jupyter_code>arsenic_list = [1.3, 1.5, 1.8, 2.6, 2.8, 3.5, 4.0, 4.8, 8, 9.5, 12, 14, 19, 23, 41, 80, 100, 110, 120, 190, 240, 250, 300, 340, 580]
matplotlib.pyplot.figure(figsize=(6,6))
#matplotlib.pyplot.bar(arsenic_list,height=2.0)
matplotlib.pyplot.boxplot(arsenic_list)
matplotlib.pyplot.show()<jupyter_output><empty_output><jupyter_text>The sample median $\hat{c}_{0.5} = 19$, which is the 13th observation ranked from smallest to largest in this sample size of 25. The binomial distribution is used to determine the 95-percent confidence interval for the true median concentration, $c_{0.5}$.<jupyter_code>print('sample median value = ',numpy.percentile(arsenic_list,50.0) )# median value using percentile function
arsenic_list = numpy.sort(arsenic_list)
index_value=13
print('\nElement',index_value,' of the sorted list is ',arsenic_list[index_value-1])<jupyter_output>sample median value =  19.0

Element 13  of the sorted list is  19.0<jupyter_text>The binomial distribution is used to determine the 95-percent confidence interval for the true median concentration, $c_{0.5}$. We obtain the critical values from the cumulative binomial distribution function by creating a list of binomial values of correct length that represents the index values, then pull from the sample actual values.<jupyter_code>n=len(arsenic_list)
p=0.5 # median is at 50%
prob = []
rank = []
for rv in range(0,n):
    x=float(rv)
    prob.append(scipy.stats.binom.cdf(x, n, p))
    rank.append(x)
    
for rv in range(0,n):
    if prob[rv]>0.05 :
        lower = arsenic_list[rv-1]
        rlo = rank[rv-1]
        break
        
for rv in range(0,n):
    if prob[rv]>0.95 :
        upper = arsenic_list[rv-1]
        rup = rank[rv-1]
        break

print('Estimate of median is ',numpy.percentile(arsenic_list,50.0))
print('Approximate 95% confidence interval estimate is the median is located between ',lower,'and', upper)

print('Actual interval probability = ',100*round(scipy.stats.binom.cdf(rup, n, p) - scipy.stats.binom.cdf(rlo, n, p),3),' %')
#print(arsenic_list[rlo-1])
#print(arsenic_list[rup+1])<jupyter_output>Estimate of median is  19.0
Approximate 95% confidence interval estimate is the median is located between  4.8 and 100.0
Actual interval probability =  92.4  %<jupyter_text>This code indicates that the lower and upper confidence intervals are at ranks 7 and 16, and that these translate to concentration values of 4.8 and 100 (the 8th and 17th values on the sorted list of concentration values. Because the sample size is relatively small (n = 25) we know that the interval will not be an exact 95-percent confidence interval. We computed the probability that the interval will contain the true value using the binomial cumulative distribution function. The result tells us that the true probability for this range is 0.924, which is close to the desired probability of 0.95. Thus, one could say the closed interval [4.8, 100] is the best approximation to a 95-percent confidence interval for the median for these observations.  Bootstrap ApproachAn alternative to the binomial distribution-based approach is to use bootstrapping. Bootstrapping (Efron and Tibshirani, 1994) is one of many computer-intensive methods that uses the observed data to represent the probability distribution from which the data were drawn, rather than assuming a normal or other theoretical distribution. Because of its robustness and coverage accuracy, the bootstrap method is commonly preferred for computing a confidence interval, especially when data are skewed or the distribution is unknown. The method also works well for data that do follow a specific distribution—in that case it returns interval endpoints very similar to those based on a distribution, for example t-intervals for the normal distribution.The bootstrap method consists of repeatedly resampling from the dataset, with replacement. The term “with replacement” means that each observation that has been selected is replaced in the dataset, which means that it can be selected multiple times in the sampling process. Two thousand to 10,000 resamples of the data are commonly used, and for many statistics it takes a small amount of computational time.<jupyter_code>import random

def one_sample(list_in,how_many):
    """Returns a list of a sample from a list of length how_many
    from list_in """
    one_sample = random.choices(list_in,k=how_many) #sampling with replacement
    return one_sample

def bootstrap_median(original_sample, pulls, replications):
    """Returns an list of bootstrapped sample medians:
    original_sample: list containing the original sample
    replications: number of bootstrap samples
    """
    #just_one_column = original_sample
    medians = [] #null list to store results
    for i in range(replications):
        this_sample = one_sample(original_sample,pulls) #get 500 pulls from one sample
        resampled_median = numpy.percentile(this_sample, 50)# get 50% percentile from the sample
        medians.append(resampled_median)
    return medians

#mylist = bootstrap_median(arsenic_list,25,200)



#matplotlib.pyplot.hist(mylist)
left_ends = []
right_ends = []
for i in range(100):
    base_sample = one_sample(arsenic_list,25)
    medians = bootstrap_median(base_sample, 25 , 100)
    left_ends.append(numpy.percentile(medians,5))
    right_ends.append(numpy.percentile(medians,95))
#matplotlib.pyplot.hist(left_ends);
#matplotlib.pyplot.hist(right_ends,alpha=0.5);
matplotlib.pyplot.figure(figsize=(8,8))

for i in range(100):
    x = []
    y = []
    x.append(left_ends[i])
    x.append(right_ends[i])
    y.append(i)
    y.append(i)
    matplotlib.pyplot.plot(x, y,color="gold")
xx = []
yy = []
xx.append(numpy.percentile(arsenic_list,50.0))
xx.append(numpy.percentile(arsenic_list,50.0))
yy.append(0)
yy.append(100)
matplotlib.pyplot.plot(xx, yy, color='red', lw=2)
matplotlib.pyplot.xlabel('Median')
matplotlib.pyplot.ylabel('Replication')
matplotlib.pyplot.title('Sample Median and Intervals of Estimates');
# count 
how_many = 0
for i in range(100):
    if left_ends[i] <= xx[0] and xx[1] <= right_ends[i] :  # count number of intervals that contain the estimate
        how_many = how_many + 1
    else:
        continue # keep going thru loop
print(how_many)
# count 
how_many = 0
for i in range(100):
    if X[i] <= xx[0] and xx[1] <= X[i+1] :  # count number of intervals that contain the estimate
        how_many = how_many + 1
#        print(X[i],X[i+1])
    else:
        continue # keep going thru loop
print(how_many)<jupyter_output>50
{'fraction_non_alphanumeric': 0.06627453828092852, 'fraction_numerical': 0.02773945247258732, 'mean_word_length': 5.083640111912826, 'pattern_counts': {'":': 0, '<': 155, '<?xml version=': 0, '>': 168, 'https://': 5, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 4, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 1, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '39', 'dir': 'jupyter-structured-clean-dedup', 'id': '32496', 'license': 'CC0-1.0', 'n_tokens_mistral': 12164, 'n_tokens_neox': 10904, 'n_words': 5302, 'path': '1-Lessons/Lesson18/.ipynb_checkpoints/ENGR-1330-Lesson18-checkpoint.ipynb', 'repo_name': 'dustykat/engr-1330-psuedo-course', 'size': '41313'}