Intervals.ipynb

{
 "metadata": {
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7-final"
  },
  "orig_nbformat": 2,
  "kernelspec": {
   "name": "python_defaultSpec_1597388694049",
   "display_name": "Python 3.7.7 64-bit (conda)"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence Interval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much of machine learning involves estimating the performance of a machine learning algorithm on unseen data.\n",
    "\n",
    "Confidence intervals are a way of quantifying the uncertainty of an estimate. They can be used to add a bounds or likelihood on a population parameter, such as a mean, estimated from a sample of independent observations from the population. Confidence intervals come from the field of estimation statistics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "    That a confidence interval is a bounds on an estimate of a population parameter.\n",
    "    That the confidence interval for the estimated skill of a classification method can be calculated directly.\n",
    "    That the confidence interval for any arbitrary population statistic can be estimated in a distribution-free way using the bootstrap.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is a Confidence Interval?\n",
    "\n",
    "A confidence interval is a bounds on the estimate of a population variable. It is an interval statistic used to quantify the uncertainty on an estimate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A confidence interval is different from a tolerance interval that describes the bounds of data sampled from the distribution. It is also different from a prediction interval that describes the bounds on a single observation. Instead, the confidence interval provides bounds on a population parameter, such as a mean, standard deviation, or similar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In applied machine learning, we may wish to use confidence intervals in the presentation of the skill of a predictive model.\n",
    "\n",
    "For example, a confidence interval could be used in presenting the skill of a classification model, which could be stated as:\n",
    "\n",
    "Given the sample, there is a 95% likelihood that the range x to y covers the true model accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The choice of 95% confidence is very common in presenting confidence intervals, although other less common values are used, such as 90% and 99.7%. In practice, you can use any value you prefer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The value of a confidence interval is its ability to quantify the uncertainty of the estimate. It provides both a lower and upper bound and a likelihood. Taken as a radius measure alone, the confidence interval is often referred to as the margin of error and may be used to graphically depict the uncertainty of an estimate on graphs through the use of error bars."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often, the larger the sample from which the estimate was drawn, the more precise the estimate and the smaller (better) the confidence interval.\n",
    "\n",
    "    Smaller Confidence Interval: A more precise estimate.\n",
    "    Larger Confidence Interval: A less precise estimate.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interval for Classification Accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Classification problems are those where a label or class outcome variable is predicted given some input data.\n",
    "\n",
    "It is common to use classification accuracy or classification error (the inverse of accuracy) to describe the skill of a classification predictive model. For example, a model that makes correct predictions of the class outcome variable 75% of the time has a classification accuracy of 75%, calculated as:\n",
    "accuracy = total correct predictions / total predictions made * 100\n",
    "1\n",
    "\t\n",
    "accuracy = total correct predictions / total predictions made * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Classification accuracy or classification error is a proportion or a ratio. It describes the proportion of correct or incorrect predictions made by the model. Each prediction is a binary decision that could be correct or incorrect. Technically, this is called a Bernoulli trial, named for Jacob Bernoulli. The proportions in a Bernoulli trial have a specific distribution called a binomial distribution. Thankfully, with large sample sizes (e.g. more than 30), we can approximate the distribution with a Gaussian."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the assumption of a Gaussian distribution of the proportion (i.e. the classification accuracy or error) to easily calculate the confidence interval.\n",
    "\n",
    "In the case of classification error, the radius of the interval can be calculated as:\n",
    "interval = z x sqrt( (error x (1 - error)) / n)\n",
    "1\n",
    "\t\n",
    "interval = z x sqrt( (error x (1 - error)) / n)\n",
    "\n",
    "In the case of classification accuracy, the radius of the interval can be calculated as:\n",
    "interval = z x sqrt( (accuracy x (1 - accuracy)) / n)\n",
    "1\n",
    "\t\n",
    "interval = z x sqrt( (accuracy x (1 - accuracy)) / n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where interval is the radius of the confidence interval, error and accuracy are classification error and classification accuracy respectively, n is the size of the sample, sqrt is the square root function, and z is the number of standard deviations from the Gaussian distribution. Technically, this is called the Binomial proportion confidence interval.\n",
    "\n",
    "Commonly used number of standard deviations from the Gaussian distribution and their corresponding significance level are as follows:\n",
    "\n",
    "    1.64 (90%)\n",
    "    1.96 (95%)\n",
    "    2.33 (98%)\n",
    "    2.58 (99%)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a model with an error of 20%, or 0.2 (error = 0.2), on a validation dataset with 50 examples (n = 50). We can calculate the 95% confidence interval (z = 1.96) as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "0.111\n"
    }
   ],
   "source": [
    "# binomial confidence interval\n",
    "from math import sqrt\n",
    "interval = 1.96 * sqrt( (0.2 * (1 - 0.2)) / 50)\n",
    "print('%.3f' % interval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then make claims such as:\n",
    "\n",
    "    The classification error of the model is 20% +/- 11%\n",
    "    The true classification error of the model is likely between 9% and 31%.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the impact that the sample size has on the precision of the estimate in terms of the radius of the confidence interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "0.078\n"
    }
   ],
   "source": [
    "# binomial confidence interval\n",
    "interval = 1.96 * sqrt( (0.2 * (1 - 0.2)) / 100)\n",
    "print('%.3f' % interval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that the confidence interval is a likelihood over a range. The true model skill may lie outside of the range."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, if we repeated this experiment over and over, each time drawing a new sample S, containing […] new examples, we would find that for approximately 95% of these experiments, the calculated interval would contain the true error. For this reason, we call this interval the 95% confidence interval estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The proportion_confint() statsmodels function an implementation of the binomial proportion confidence interval.\n",
    "\n",
    "By default, it makes the Gaussian assumption for the Binomial distribution, although other more sophisticated variations on the calculation are supported. The function takes the count of successes (or failures), the total number of trials, and the significance level as arguments and returns the lower and upper bound of the confidence interval.\n",
    "\n",
    "The example below demonstrates this function in a hypothetical case where a model made 88 correct predictions out of a dataset with 100 instances and we are interested in the 95% confidence interval (provided to the function as a significance of 0.05)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "lower=0.816, upper=0.944\n"
    }
   ],
   "source": [
    "from statsmodels.stats.proportion import proportion_confint\n",
    "lower, upper = proportion_confint(88, 100, 0.05)\n",
    "print('lower=%.3f, upper=%.3f' % (lower, upper))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nonparametric Confidence Interval\n",
    "\n",
    "Often we do not know the distribution for a chosen performance measure. Alternately, we may not know the analytical way to calculate a confidence interval for a skill score."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In these cases, the bootstrap resampling method can be used as a nonparametric method for calculating confidence intervals, nominally called bootstrap confidence intervals.\n",
    "\n",
    "The bootstrap is a simulated Monte Carlo method where samples are drawn from a fixed finite dataset with replacement and a parameter is estimated on each sample. This procedure leads to a robust estimate of the true population parameter via sampling.\n",
    "\n",
    "We can demonstrate this with the following pseudocode."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "statistics = []\n",
    "for i in bootstraps:\n",
    "\tsample = select_sample_with_replacement(data)\n",
    "\tstat = calculate_statistic(sample)\n",
    "\tstatistics.append(stat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The procedure can be used to estimate the skill of a predictive model by fitting the model on each sample and evaluating the skill of the model on those samples not included in the sample. The mean or median skill of the model can then be presented as an estimate of the model skill when evaluated on unseen data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Confidence intervals can be added to this estimate by selecting observations from the sample of skill scores at specific percentiles.\n",
    "\n",
    "Recall that a percentile is an observation value drawn from the sorted sample where a percentage of the observations in the sample fall. For example, the 70th percentile of a sample indicates that 70% of the samples fall below that value. The 50th percentile is the median or middle of the distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we must choose a significance level for the confidence level, such as 95%, represented as 5.0% (e.g. 100 – 95). Because the confidence interval is symmetric around the median, we must choose observations at the 2.5th percentile and the 97.5th percentiles to give the full range.\n",
    "\n",
    "We can make the calculation of the bootstrap confidence interval concrete with a worked example.\n",
    "\n",
    "Let’s assume we have a dataset of 1,000 observations of values between 0.5 and 1.0 drawn from a uniform distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bootstrap confidence intervals\n",
    "import numpy as np\n",
    "from numpy.random import seed\n",
    "from numpy.random import rand\n",
    "from numpy.random import randint\n",
    "from numpy import mean\n",
    "from numpy import median\n",
    "from numpy import percentile\n",
    "# generate dataset\n",
    "dataset = 0.5 + rand(1000) * 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "0.7509369652161508"
     },
     "metadata": {},
     "execution_count": 21
    }
   ],
   "source": [
    "np.mean(dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will perform the bootstrap procedure 100 times and draw samples of 1,000 observations from the dataset with replacement. We will estimate the mean of the population as the statistic we will calculate on the bootstrap samples. This could just as easily be a model evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bootstrap\n",
    "scores = list()\n",
    "for _ in range(100):\n",
    "\t# bootstrap sample\n",
    "\tindices = randint(0, 1000, 1000)\n",
    "\tsample = dataset[indices]\n",
    "\t# calculate and store statistic\n",
    "\tstatistic = mean(sample)\n",
    "\tscores.append(statistic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(625, 1000, 100, 0.7432532429545529)"
     },
     "metadata": {},
     "execution_count": 22
    }
   ],
   "source": [
    "len(np.unique(indices)), len(indices), len(scores), scores[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have a sample of bootstrap statistics, we can calculate the central tendency. We will use the median or 50th percentile as we do not assume any distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "median=0.751\n"
    }
   ],
   "source": [
    "print('median=%.3f' % median(scores))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then calculate the confidence interval as the middle 95% of observed statistical values centered around the median."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate 95% confidence intervals (100 - alpha)\n",
    "alpha = 5.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, the desired lower percentile is calculated based on the chosen confidence interval. Then the observation at this percentile is retrieved from the sample of bootstrap statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate lower percentile (e.g. 2.5)\n",
    "lower_p = alpha / 2.0\n",
    "# retrieve observation at lower percentile\n",
    "lower = max(0.0, percentile(scores, lower_p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do the same thing for the upper boundary of the confidence interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate upper percentile (e.g. 97.5)\n",
    "upper_p = (100 - alpha) + (alpha / 2.0)\n",
    "# retrieve observation at upper percentile\n",
    "upper = min(1.0, percentile(scores, upper_p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "50th percentile (median) = 0.751\n2.5th percentile = 0.742\n97.5th percentile = 0.757\n"
    }
   ],
   "source": [
    "print('50th percentile (median) = %.3f' % median(scores))\n",
    "print('%.1fth percentile = %.3f' % (lower_p, lower))\n",
    "print('%.1fth percentile = %.3f' % (upper_p, upper))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running the example summarizes the distribution of bootstrap sample statistics including the 2.5th, 50th (median) and 97.5th percentile.\n",
    "\n",
    "We can then use these observations to make a claim about the sample distribution, such as:\n",
    "\n",
    "There is a 95% likelihood that the range 0.742 to 0.757 covers the true statistic mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Example for ML Model\n",
    "import numpy\n",
    "from pandas import read_csv\n",
    "from sklearn.utils import resample\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.metrics import accuracy_score\n",
    "from matplotlib import pyplot\n",
    "# load dataset\n",
    "data = read_csv('pima-indians-diabetes.data.csv', header=None)\n",
    "values = data.values\n",
    "# configure bootstrap\n",
    "n_iterations = 1000\n",
    "n_size = int(len(data) * 0.50)\n",
    "# run bootstrap\n",
    "stats = list()\n",
    "for i in range(n_iterations):\n",
    "\t# prepare train and test sets\n",
    "\ttrain = resample(values, n_samples=n_size)\n",
    "\ttest = numpy.array([x for x in values if x.tolist() not in train.tolist()])\n",
    "\t# fit model\n",
    "\tmodel = DecisionTreeClassifier()\n",
    "\tmodel.fit(train[:,:-1], train[:,-1])\n",
    "\t# evaluate model\n",
    "\tpredictions = model.predict(test[:,:-1])\n",
    "\tscore = accuracy_score(test[:,-1], predictions)\n",
    "\tprint(score)\n",
    "\tstats.append(score)\n",
    "# plot scores\n",
    "pyplot.hist(stats)\n",
    "pyplot.show()\n",
    "# confidence intervals\n",
    "alpha = 0.95\n",
    "p = ((1.0-alpha)/2.0) * 100\n",
    "lower = max(0.0, numpy.percentile(stats, p))\n",
    "p = (alpha+((1.0-alpha)/2.0)) * 100\n",
    "upper = min(1.0, numpy.percentile(stats, p))\n",
    "print('%.1f confidence interval %.1f%% and %.1f%%' % (alpha*100, lower*100, upper*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prediction Interval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://machinelearningmastery.com/time-series-forecast-uncertainty-using-confidence-intervals-python/\n",
    "\n",
    "Time series forecast models can both make predictions and provide a prediction interval for those predictions.\n",
    "\n",
    "Prediction intervals provide an upper and lower expectation for the real observation. These can be useful for assessing the range of real possible outcomes for a prediction and for better understanding the skill of the model\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tolerance Interval\n",
    "\n",
    "https://machinelearningmastery.com/statistical-tolerance-intervals-in-machine-learning/\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be useful to have an upper and lower limit on data.\n",
    "\n",
    "These bounds can be used to help identify anomalies and set expectations for what to expect. A bound on observations from a population is called a tolerance interval. A tolerance interval comes from the field of estimation statistics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A tolerance interval is different from a prediction interval that quantifies the uncertainty for a single predicted value. It is also different from a confidence interval that quantifies the uncertainty of a population parameter such as a mean. Instead, a tolerance interval covers a proportion of the population distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After completing this tutorial, you will know:\n",
    "\n",
    "    That statistical tolerance intervals provide a bounds on observations from a population.\n",
    "    That a tolerance interval requires that both a coverage proportion and confidence be specified.\n",
    "    That the tolerance interval for a data sample with a Gaussian distribution can be easily calculated.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bounds on Data\n",
    "\n",
    "It is useful to put bounds on data.\n",
    "\n",
    "For example, if you have a sample of data from a domain, knowing the upper and lower bound for normal values can be helpful for identifying anomalies or outliers in the data.\n",
    "\n",
    "For a process or model that is making predictions, it can be helpful to know the expected range that sensible predictions may take.\n",
    "\n",
    "Knowing the common range of values can help with setting expectations and detecting anomalies.\n",
    "\n",
    "The range of common values for data is called a tolerance interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What Are Statistical Tolerance Intervals?\n",
    "\n",
    "The tolerance interval is a bound on an estimate of the proportion of data in a population.\n",
    "\n",
    "The interval is limited by the sampling error and by the variance of the population distribution. Given the law of large numbers, as the sample size is increased, the probabilities will better match the underlying population distribution.\n",
    "\n",
    "Below is an example of a stated tolerance interval:\n",
    "\n",
    "The range from x to y covers 95% of the data with a confidence of 99%.\n",
    "\n",
    "If the data is Gaussian, the interval can be expressed in the context of the mean value; for example:\n",
    "\n",
    "x +/- y covers 95% of the data with a confidence of 99%.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A tolerance interval is defined in terms of two quantities:\n",
    "\n",
    "    Coverage: The proportion of the population covered by the interval.\n",
    "    Confidence: The probabilistic confidence that the interval covers the proportion of the population.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How to Calculate Tolerance Intervals\n",
    "\n",
    "The size of a tolerance interval is proportional to the size of the data sample from the population and the variance of the population.\n",
    "\n",
    "There are two main methods for calculating tolerance intervals depending on the distribution of data: parametric and nonparametric methods.\n",
    "\n",
    "    Parametric Tolerance Interval: Use knowledge of the population distribution in specifying both the coverage and confidence. Often used to refer to a Gaussian distribution.\n",
    "    Nonparametric Tolerance Interval: Use rank statistics to estimate the coverage and confidence, often resulting less precision (wider intervals) given the lack of information about the distribution.\n",
    "\n",
    "Tolerance intervals are relatively straightforward to calculate for a sample of independent observations drawn from a Gaussian distribution. We will demonstrate this calculation in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tolerance Interval for Gaussian Distribution\n",
    "\n",
    "In this section, we will work through an example of calculating the tolerance intervals on a data sample.\n",
    "\n",
    "First, let’s define our data sample. We will create a sample of 100 observations drawn from a Gaussian distribution with a mean of 50 and a standard deviation of 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate dataset\n",
    "from numpy.random import randn\n",
    "data = 5 * randn(100) + 50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During the example, we will assume that we are unaware of the true population mean and standard deviation, and that these values must be estimated.\n",
    "\n",
    "Because the population parameters have to be estimated, there is additional uncertainty. For example, for a 95% coverage, we could use 1.96 (or 2) standard deviations from the estimated mean as the tolerance interval. We must estimate the mean and standard deviation from the sample and take this uncertainty into account, therefore the calculation of the interval is slightly more complex.\n",
    "\n",
    "Next, we must specify the number of degrees of freedom. This will be used in the calculation of critical values and in the calculation of the interval. Specifically, it is used in the calculation of the standard deviation.\n",
    "\n",
    "Remember that the degrees of freedom are the number of values in the calculation that can vary. Here, we have 100 observations, therefore 100 degrees of freedom. We do not know the standard deviation, therefore it must be estimated using the mean. This means our degrees of freedom will be (N – 1) or 99."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# specify degrees of freedom\n",
    "n = len(data)\n",
    "dof = n - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we must specify the proportional coverage of the data. In this example, we are interested in the middle 95% of the data. The proportion is 95. We must shift this proportion so that it covers the middle 95%, that is from 2.5th percentile to the 97.5th percentile.\n",
    "\n",
    "We know that the critical value for 95% is 1.96 given that we use it so often; nevertheless, we can calculate it directly in Python given the percentage 2.5% of the inverse survival function. This can be calculated using the norm.isf() SciPy function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import chi2\n",
    "from scipy.stats import norm\n",
    "# specify data coverage\n",
    "prop = 0.95\n",
    "prop_inv = (1.0 - prop) / 2.0\n",
    "gauss_critical = norm.isf(prop_inv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to calculate the confidence of the coverage. We can do this by retrieving the critical value from the Chi Squared distribution for the given number of degrees of freedom and desired probability. We can use the chi2.isf() SciPy function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# specify confidence\n",
    "prob = 0.99\n",
    "chi_critical = chi2.isf(q=prob, df=dof)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have all of the pieces to calculate the Gaussian tolerance interval. The calculation is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "interval = np.sqrt((dof * (1 + (1/n)) * gauss_critical^2) / chi_critical)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where dof is the number of degrees of freedom, n is the size of the data sample, gauss_critical is the critical value, such as 1.96 for 95% coverage of the population, and chi_critical is the Chi Squared critical value for the desired confidence and degrees of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interval = sqrt((dof * (1 + (1/n)) * gauss_critical**2) / chi_critical)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "Gaussian critical value: 1.960 (coverage=95%)\nChi-Squared critical value: 69.230 (prob=99%, dof=99)\nTolerance Interval: 2.355\n47.95 to 52.66 covers 95% of data with a confidence of 99%\n"
    }
   ],
   "source": [
    "# parametric tolerance interval\n",
    "from numpy.random import seed\n",
    "from numpy.random import randn\n",
    "from numpy import mean\n",
    "from numpy import sqrt\n",
    "from scipy.stats import chi2\n",
    "from scipy.stats import norm\n",
    "# seed the random number generator\n",
    "seed(1)\n",
    "# generate dataset\n",
    "data = 5 * randn(100) + 50\n",
    "# specify degrees of freedom\n",
    "n = len(data)\n",
    "dof = n - 1\n",
    "# specify data coverage\n",
    "prop = 0.95\n",
    "prop_inv = (1.0 - prop) / 2.0\n",
    "gauss_critical = norm.isf(prop_inv)\n",
    "print('Gaussian critical value: %.3f (coverage=%d%%)' % (gauss_critical, prop*100))\n",
    "# specify confidence\n",
    "prob = 0.99\n",
    "chi_critical = chi2.isf(q=prob, df=dof)\n",
    "print('Chi-Squared critical value: %.3f (prob=%d%%, dof=%d)' % (chi_critical, prob*100, dof))\n",
    "# tolerance\n",
    "interval = sqrt((dof * (1 + (1/n)) * gauss_critical**2) / chi_critical)\n",
    "print('Tolerance Interval: %.3f' % interval)\n",
    "# summarize\n",
    "data_mean = mean(data)\n",
    "lower, upper = data_mean-interval, data_mean+interval\n",
    "print('%.2f to %.2f covers %d%% of data with a confidence of %d%%' % (lower, upper, prop*100, prob*100))"
   ]
  }
 ]
}