Comprehensive Confidence Intervals for Python Developers

Confidence interval is uncertainty in summary statistic represented as a range. In the other words, it is a range of values we are fairly sure our true value lies in. For example: I am 95% confident that the population mean falls between 8.76 and 15.88 $\rightarrow$ (12.32 $\pm$ 3.56)

Confidence interval tells you how confident you can be that the results from a poll or survey reflect what you would expect to find if it were possible to survey the entire population. It is difficult to obtain measurement data of an entire data set (population) due to limited resource & time. Your best shot is to survey a small fraction (samples) of the entire data set, and pray that your sample data represents the population reasonably well.

Sample data may not be a good representation of a population by numerous factors (Ex: bias), and as a result, uncertainty is always introduced in any estimations derived from sample data. Due to the uncertainty involved with sample data, any statistical estimation needs to be delivered in a range, not in a point estimate.

How well a sample statistic estimates an underlying population parameter is always an issue (Population vs. Samples). A confidence interval addresses this issue by providing a range of values, which is likely to contain the population parameter of interest within the range of uncertainty.

1. Understanding confidence interval with analogy¶

If you've taken a science class with lab reports in your highschool or college, you probably had to include measurement error in your lab reports. For example, if you were asked to measure the length of a paper clip with a ruler, you have to include $\pm0.5 \,\text{cm}$ or $\pm0.05\,\text{cm}$ (depending on the spacing of tick marks) to account for the measurement error that shows the precision of your measuring tool.

Based on figure (1), the paper clip seems to be about 2.7 cm long, but we don't know for sure because the tickmarks in the ruler is not precise enough to measure decimal length. However, I can tell with 100% confidence that the paper clip has a length between 2 ~ 3 cm, because the clip is between the 2 cm and 3 cm tickmarks. You record the length of the paper clip in a range, instead of a point estimate, to account for the uncertainty introduced by the limitation of the measuring tool.

Figure 1: Measurement error in ruler

Similar idea can be applied to a confidence interval of mean. You want to obtain a mean of a whole data set (population), but you can measure values of only a small fraction (samples) of the whole data set. This boils down to the traditional issue of Population vs Samples, due to the cost of obtaining measurement data of a large data set. Uncertainty is introduced in your samples, because you don't know if your samples are 100% representative of the population, free of bias. Therefore, you deliver your conclusion in a range, not in a point estimate, to account for the uncertainty.

Example 1: Uncertainty in rock porosity

(Borrowed from Dr. Michael Pyrcz's Geostatistics class)

A reservoir engineer in the oil & gas industry wants to know the rock porosity of a formation to estimate the total oil reserve 9,500 ft underground. Due to the high cost of obtaining rock core samples from the deep formations, he could acquire only 12 rock core samples. Since the uncertainty of a point estimation scales inversely with a sample size, his estimation is subject to non-negligible uncertainty. He obtains 14.5% average rock porosity with 4.3% standard deviation. Executives in the company wants to know the worst-case scenario (P10) and the best-case scenario (P90) to make business decisions. You can convey your estimation of average porosity with uncertainty by constructing the confidence interval of mean.

Assuming that you have a reason to believe that the rock porosity follows normal distribution, you can construct its 80% confidence interval, with the procedure described below:

The above range of uncertainty was acquired from the 12 rock core samples. In the worst-case scenario, the rock formation at 9,500 ft underground has 12.8% porosity. In the best-case scenario, the oil reservoir has 16.2% porosity. The same procedures can be applied for the core samples collected at different depths, which give us the confidence interval plot of rock porosities shown in figure (2).

Figure 2: Confidence interval of core samples porosities along depths

Example 2: Purity of methamphetamine (crystal) in Breaking Bad

21 batches of crystal cooked by Mr. White shows 99.1% average purity with 3% standard deviation. 18 batches of crystal cooked by Mr. Pinkman shows 96.2% average purity with 4% standard deviation. Does Mr. White always cook better crystal than Mr. Pinkman, or is it possible for Mr. Pinkman to beat Mr. White in purity of cooked crystals, by luck? We can construct 95% confidence interval assuming normal distribution, with the procedure described below:

There's a small overlap between the confidence intervals of Mr. White's and Mr. Pinkman's. Although it is true that Mr. White is a better cooker, Mr. Pinkman can cook a purer batch of crystals by a small chance, if he has the luck. Comparing the means of two sample data sets is closely related to constructing confidence interval of difference in mean.

Figure 3: Overlap in the 95% confidence interval of two samples

Example 3: Uncertainty in oil production forecast

A production engineer in the oil & gas industry wants to know the worst-case scenario (P10) and the best-case scenario (P90) of hydrocarbon production forecast. In petroleum engineering, we use a technique called Decline Curve Analysis (DCA) to project future hydrocarbon production. It is important to quantify the uncertanties of your DCA model, as the uncertainty in oil production can be as large as millions of dollars worth; executives in a company make business decisions based on P10 and P90 values. There are mainly three models for DCA: exponential, hyperbolic, and harmonic. For demonstration purpose, hyperbolic model will be used here. Hyperbolic decline curve can be defined as:

It is a non-linear regression problem with three parameters to optimize: $Di$, $q_i$, and $b$. The independent variable is time $t$. We can construct the confidence interval of regression model to get P10 and P90 values. Since the analytical solution for the confidence interval of non-linear model is complicated, and the data is not normally distributed, we use non-parametric numerical alternative, bootstrap, to obtain the following uncertainty forecast model:

Figure 4: Uncertainty in decline curve analysis forecast

2. Key takeaways¶

3. Population vs. samples¶

Confidence interval describes the amount of uncertainty associated with a sample estimate of a population parameter. One needs to have a good understanding of the difference between samples and population to understand the necessity of delivering statistical estimations in a range, a.k.a. confidence interval.

Figure 5: Population vs samples

Population: data set that contains all members of a specified group. Ex: ALL people living in the US.

Samples: data set that contains a part, or a subset, of a population Ex: SOME people living in the US.

Let's say that you are conducting a phone-call survey to investigate the society's perception of The Affordable Care Act (“Obamacare”). Since you can't call all 327.2 million people (population) in the US, you call about 1,000 people (samples). Your poll showed that 59% of the registered voters support Obamacare. This does not agree with the actual survey conducted in 2018; 53% favorable, 42% unfavorable (source). What could be the source of error?

Since (formal) president Obama is a member of the Democratic Party, the voters' response can be affected by their political preference. How could you tell that the 1,000 people you called happened to be mostly Democrats, who's more likely to support Obama's policy, because they share similar political view? The samples you collected could have been biased, but you don't that know for sure. Of course, the voters' response could be affected by many other factors like race, age, place of residence, or financial status. The idea is that, there will always be uncertainty involved with your estimation, because you don't have an access to the entire population.

Confidence interval is a technique that quantifies the uncertainty when estimating a population parameter from samples.

Pythonic Tip: Difference between Numpy variance and Pandas variance

Different libraries make different assumption about an input array. The default value of ddof is different for Pandas and Numpy, resulting in different variance. ddof represent degrees of freedom, and setting ddof=True or ddof=1 tells the variance function to calculate sample variance by accounting for the bias factor $n-1$ (recall that in Python, True==1.) Remember that there is a distinction between Population variance ($\sigma^2$) vs. Sample variance ($s^2$). If you are confused which library is computing which variance (sample or population), just remember this: whatever library you are using, use ddof=True or ddof=1 to compute sample variance, and use ddof=False or ddof=0 to compute population variance.

4. Confidence interval of normal distribution¶

Computing confidence interval of a statistic depends on two factors: type of statistic, and type of sample distribution. As explained above, different formulas exist for different type of statistics (Ex: mean, std, variance), and different methods (Ex: boostrapping, credible interval, Box-Cox transformation) are used for non-normal data set.

We will cover confidence interval of mean, difference in mean and variance.

4.1. Confidence interval of mean¶

Confidence interval of mean is used to estimate the population mean from sample data and quantify the related uncertainty. Consider the following figure:

Figure 6: Distribution of population and C.I. of mean

In figure (6), assume that the population is normally distributed. Since we don't have an access to the entire population, we have to guess the population mean (unknown) to the best of our ability using sample data set. We do this by computing the sample mean and constructing its 95% confidence interval. Note that the popular choices of confidence level are: 90%, 95%, and 99%

Assuming normality of population, its sample means are also normally distributed. Let's say that you have a population, and you draw small fractions of it $N$ times. Then, the computed means of $N$ sample sets $\boldsymbol{\mu}=(\mu_1, \mu_2,..., \mu_{N-1}, \mu_N)$ is normally distributed as shown in figure (7). Their confidence intervals are represented as the black horizontal arrows:.

Figure 7: Distribution of sample mean and its C.I.

You can see that the confidence interval of $\mu_5$ does NOT include the green vertical dashed line, 12.31. Let's assume that 12.31 is the true population mean (we never know if this is the actual population mean or not, but let's assume). If we get $\mu_5$ and its confidence interval as our estimation of the population mean, then our estimation is wrong. There is a 5% chance of this happening, because we set our confidence level as 95%. Note that the width of the confidence intervals (black horizontal arrows) depend on the sample size, as shown in eq (1)

The grey area of figure (6) is essentially equivalent to the grey area of figure (7). $\mu_1$ = 12.32 is the sample mean, and $\pm$ 3.56 is the uncertainty related to the sample mean with 95% confidence. The uncertainty is a product of distribution score and standard error of mean. Distribution score essentially tells how many standard error are the limits (8.76 and 15.88) away from the center (12.32). Choosing larger confidence level results in larger confidence interval. This increases the grey area in figure (6) and figure (7).

We convey 95% confidence interval of mean like this:

I am 95% confident that the population mean falls between 8.76 and 15.88. If I sample data 20 times, 19 times the sample mean will fall between 8.76 ~ 15.88, but expect that I will be wrong 1 time.

Notes: Distribution of various statistics

Different statistics exhibit different distributions. Normality of samples does not guarantee normality of its statistics. When the samples are normally distributed, their means are normally distributed, but their variances are chi-square $\chi^2$ distributed. More discussion about the distribution of variance and $\chi^2$ distribution is covered below. Note that these assumptions are invalid when samples are non-normal.

Source Code For The Figure

                
                    from scipy import stats
                    import matplotlib.pyplot as plt
                    import numpy as np

                    df_values = [1, 2, 6, 9]
                    linestyles = ['-', '--', ':', '-.']
                    normal_params = [(10, 1), (11, 1), (10, 2), (10, 3)]
                    x = np.linspace(-1, 20, 1000)

                    plt.style.use('seaborn-whitegrid')
                    fig, ax = plt.subplots(1, 2, figsize=(13.3, 5))
                    fig.tight_layout()
                    plt.subplots_adjust(left=0.09, right=0.96, bottom=0.12, top=0.93)

                    for df, norm_p, ls in zip(df_values, normal_params, linestyles):
                        ax[1].plot(x, stats.chi2.pdf(x, df, loc=0, scale=1),
                                   ls=ls, c='black', label=r'Degrees of freedom$=%i$' % df)
                        ax[0].plot(x, stats.norm.pdf(x, loc=norm_p[0], scale=norm_p[1]), 
                                   ls=ls, c='black', label='Mean = %d, ' % norm_p[0] + 'Std = %s' % norm_p[1])

                    ax[0].set_xlim(4, 16)
                    ax[0].set_ylim(-0.025, 0.525)
                    ax[0].set_xlabel('$x$', fontsize=20)
                    ax[0].set_ylabel(r'Probability', fontsize=20)
                    ax[0].set_title(r'Distribution of means: normal distribution', fontsize=20)
                    ax[0].legend(loc='upper left', fontsize=16, framealpha=1, frameon=True)

                    ax[1].set_xlim(0, 10)
                    ax[1].set_ylim(-0.025, 0.525)
                    ax[1].set_xlabel('$\chi^2$', fontsize=20)
                    ax[1].set_title(r'Distribution of variances: $\chi^2$ distribution', fontsize=20)
                    ax[1].legend(loc='best', fontsize=16, framealpha=1, frameon=True)
                
            

If sample data is normal or normal-like distributed, we almost always assume t-distribution to compute confidence interval, as explained below. Then, the confidence interval of mean has the following analytical solution:

Recall that when computing $s$, correction factor ($n-1$) is applied to account for sample bias, as explained above. Pay close attention to the standard error $\frac{s}{\sqrt{(n)}}$. As the sample size $n$ increases, the standard error decreases, reducing the range of confidence interval. This is intuitive in a sense that, the more samples we have, the less uncertainty we have with our statistical estimation. The length of the black horizontal arrows in figure (7) depends on the sample size. The larger the sample size, the narrower the width of arrows, and vice versa.

Pythonic Tip: Computing confidence interval of mean with SciPy

We can compute confidence interval of mean directly from using eq (1). Recall to pass ddof=1 to make sure to compute sample standard deviation $s$, not population standard deviation $\sigma$, as explained above. We will draw random samples from normal distribution using np.random.normal(). Note that loc is for population mean, and scale is for population standard deviation, and size is for number of samples to draw.

Or we can compute with scipy.stats.t.interval(). Note that you don't divide alpha by 2, because the function does that for you. Also note that the standard error of mean $\frac{s}{\sqrt{n}}$ can be computed with scipy.stats.sem()

Note the default value of loc=0 and scale=1. This will assume sample mean $\mu$ to be 0, and standard error $\frac{s}{\sqrt{n}}$ to be 1, which assumes standard normal distribution of mean = 0 and standard deviation = 1. This is NOT what we want.

4.2. Confidence interval of difference in mean¶

Confidence interval of difference in mean is not very useful by itself. But it is important to understand how it works, because it forms the basis of one of the most widely used hypothesis test: t-test.

Often we are interested in knowing if two distributions are significantly different. In the other words, we want to know if two sample data sets came from the same population by comparing central tendency of populations. A standard approach is to check if the sample means are different. However, this is a misleading approach in a sense that the means of samples are almost always different, even if the difference is microscopic. More useful would be to estimate the difference in a range to account for uncertainty, and compute probability that it is big enough to be of practical importance. T-test checks if the difference is "close enough" to zero by computing the confidence interval of difference in means.

T-test hypothesis

Note that the above hypothesis tests whether the mean of one group is significantly DIFFERENT from the mean of the other group; we are using two-tailed test. This does not check if the mean of one group is significantly GREATER than the mean of the other group, which uses one-tailed test.

Figure 8: Distributions of samples

In figure (8), $\mu$ represents the sample mean. If two sample data sets are from the same population, the distribution of means will be similar "enough". If not, they will be "significantly" different. It can be visually inspected by the area of overlap. The larger the overlap, the bigger the chance of the two distributions originating from the same population.

The more robust way to compare sample means would be to construct the confidence interval of difference in means. If the two samples came from the same population, they should have the similar "enough" means. Their difference should be close to zero and satisfy (or fail to reject) the null hypothesis $H_0: \mu_1 - \mu_2 = 0$ within a range of uncertainty. Consider the following figure:

Figure 9: Distribution of difference in means

In figure (9), the calculated difference in sample means is $\mu_1 - \mu_2 = 1.00$. We deliver the uncertainty related to our estimation of difference in sample means by constructing its 95% confidence interval $[$-1.31 ~ 3.31$]$. Since the null hypothesis $H_0: \mu_1 - \mu_2 = 0$ is within the 95% confidence interval (grey shaded area), we accept the null hypothesis; we conclude that the samples have the same means within the uncertainty.

However, if the null hypothesis is not within the confidence interval and falls in the 2.5% outliers zone, we reject the null hypothesis and accept the alternate hypothesis $H_1: \mu_1 - \mu_2 \neq 0$. In the other words, we conclude that the sample means are significantly different.

Three variations of confidence interval of difference in means

There are three variations of t-test, and therefore there are three variations of confidence interval of difference in means. The difference & application of the three variations are really well-explained in Wikipedia (one of the few that are actually easy to understand, with minimum jargons.)

1. Independent (unpaired) samples, equal variance - Student's t-interval

2. Independent (unpaired) samples, unequal variance - Welch's t-interval

3. Dependent (paired) samples

Recall that all t-tests assume normality of data. However, they are pretty robust to non-normality as long as the deviation from normality isn't large. Visualize your distributions to test this. Robustness of t-test to non-normality is discussed in detail below.

4.2.1. Independent (unpaired) samples, equal variance - student's t-interval¶

When you have a reason to believe that samples have nearly equal variances, you can use student's t-test to check if difference in means are significantly different. Note that student's t-test works pretty well even with unequal variances as long as sample sample sizes are equal or nearly equal, and sample sizes are not tiny.

However, it is recommended to always use Welch's t-test by assuming unequal variances, as explained below. Use student's t-test if you are ABSOLUTELY sure that the population variances are nearly equal.

Confidence interval of difference in mean assuming equal variance (student's t-interval) can be calculated as follows:

The formula for the pooled standard deviation $s_p$ looks a bit overwhelming, but its just an weighted average standard deviation of two samples, with bias correction factor $n_i-1$ for each sample. Recall that student's t-test assumes equal variances of two samples. You calculate what is assumed to be the common variance (=pooled variance, $s_p^2$) by computing the weighted average from each sample's variance.

In eq (4), $t$-score depends on significance level $\alpha$ and degrees of freedom $df$. In student's t-test, which assumes equal variance:

Pythonic Tip: Computing student's t-interval

Unfortunately, SciPy doesn't support computing confidence intereval of difference in mean separately. It is incorporated into computing t-statistic and p-value of t-test, but users can't access its underlying confidence interval. Note that in R, users have access to the CI of difference in means. We can compute CI of difference in means assuming equal variance with eq (4). Don't forget to compute sample variance, instead of population variance by setting ddof=1 as explained above.

The 95% confidence interval of difference in means has 0 within its interval. This means that the null hypothesis, $H_0: \mu_1 - \mu_2 = 0$ in figure (9), falls within the interval and we fail to reject the null hypothesis. We conclude that the sample means are not significantly different.

We can confirm this by running a formal hypothesis testing with scipy.stats.ttest_ind(), and setting equal_var=True. Note that this assumes independent t-test with pooled variance, which is equivalent to student's t-test.

The computed pvalue=0.229 is bigger than the significance level of alpha = 0.05, and therefore we fail to reject the null hypothesis, which is consistent with the conclusion drawn from the confidence interval of difference in mean.

Checking results with R:

    
        a <- c(12.9, 10.2, 7.4, 7.0, 10.5, 11.9, 7.1, 9.9, 14.4, 11.3)

        b <- c(10.2, 6.9, 10.9, 11.0, 10.1, 5.3, 7.5, 10.3, 9.2, 8.8)

        t.test(a, b, var.equal = TRUE)

        #   Two Sample t-test

        # data:  a and b
        # t = 1.2453, df = 18, p-value = 0.229
        # 95 percent confidence interval:
        #  -0.8520327  3.3320327
        # sample estimates:
        # mean of x mean of y 
        #     10.26      9.02 
    

4.2.2. Independent (unpaired) samples, unequal variance - Welch's t-interval¶

When comparing central tendency of normal distributions, it is safer, and therefore recommended to always use Welch's t-test, which assumes unequal variances of samples, as explained below. Equal variance t-test is not robust when population variances are different, but unequal variances are robust even when population variances are equal.

Confidence interval of difference in mean assuming unequal variance (Welch's t-interval) can be calculated as follows:

The formula is very similar to student's t-interval. There are two main differences:

1. We use each sample's own variance $s_1^2$ and $s_2^2$, instead of pooled (weighted average) variance $s_p^2$.

2. Degrees of freedom $df$ is computed with eq (7).

Pythonic Tip: Computing Welch's t-interval

The procedure is very similar to Computing student's t-interval. We will compute confidence interval of difference in mean assuming unequal variance, with eq (6). Although Scipy supports computing t-statistic for Welch's t-test, it doesn't support a function that allows us to compute Welch's t-interval. We will have to write our own codes to compute it. Don't forget to compute sample variance, instead of population variance by setting ddof=1 as explained above.

The 95% confidence interval of difference in means has 0 within its interval. This means that the null hypothesis, $H_0: \mu_1 - \mu_2 = 0$ in figure (9), falls within the interval and we fail to reject the null hypothesis. We conclude that the sample means are not significantly different.

We can confirm this by running a formal hypothesis testing with scipy.stats.ttest_ind(), and setting equal_var=False. Note that this assumes independent t-test with pooled variance, which is equivalent to student's t-test.

The computed pvalue=0.230 is bigger than the significance level of alpha = 0.05, and therefore we fail to reject the null hypothesis, which is consistent with the conclusion drawn from the confidence interval of difference in mean.

Checking results with R:

    
        a <- c(12.9, 10.2, 7.4, 7.0, 10.5, 11.9, 7.1, 9.9, 14.4, 11.3)

        b <- c(10.2, 6.9, 10.9, 11.0, 10.1, 5.3, 7.5, 10.3, 9.2, 8.8)

        t.test(a, b, var.equal = FALSE)

        #   Welch Two Sample t-test

        # data:  a and b
        # t = 1.2453, df = 16.74, p-value = 0.2302
        # alternative hypothesis: true difference in means is not equal to 0
        # 95 percent confidence interval:
        #  -0.8633815  3.3433815
        # sample estimates:
        # mean of x mean of y 
        #     10.26      9.02 

    

4.2.3. Dependent (paired) samples - Paired t-interval¶

This test is used when the samples are dependent; that is, when there is only one sample that has been tested twice (repeated measures) or when there are two samples that have been matched or "paired" (paired or unpaired? read below.)

Confidence interval of difference in means assuming paired samples can be calculated as follows:

The equation is very similar to eq (1), except that we are computing mean and standard deviation of differences between before & after state of test subjects. Let's try to understand this with an example.

A school develops a tutoring program to improve the SAT scores of high school students. A school requires students to take tests before & after tutoring, and checks if the tutoring had a significant impact on the SAT scores of students. Because the test subjects are compared to themselves, not anyone elses, the measurements taken before & after the training are not independent.

To compute dependent t-interval, we compute differences of test scores before & after tutoring:

We find $\bar{d}$ = -193.0, and $s_d$ = 161.7. These values are plugged into eq (8). Degrees of freedom $df$ for dependent t-interval can be computed with:

Unlike independent t-test, in which two samples can have different sample sizes $n_1$ and $n_2$, depedent t-test has only one sample size, because the test subjects are compared to themselves.

Also note that dependent t-test assumes difference of test scores to be normally distributed, not test scores of students themselves. But as long as the test scores are normally distributed, the difference of test scores will also be normally distributed due to the property of normal distributions.

Pythonic Tip: Computing paired t-interval

Although Scipy supports computing t-statistic for dependent t-test, it doesn't support a function that allows us to compute dependent t-interval. We will have to write our own codes to compute it. You can compute it with eq (8) Don't forget to compute sample standard devaition, instead of population standard deviation by setting ddof=1 as explained above.

The 95% confidence interval of difference in means for dependent samples does not have 0 within its interval. This means that the null hypothesis, $H_0: \mu_1 - \mu_2 = 0$ in figure (9), does not fall within the interval. Instead, our estimation falls within the 2.5% outlier zone on the left, $H_1: \mu_1 - \mu_2 \neq 0$. We reject the null hypothesis $H_0$, and accept the alternate hypothesis $H_1$. We conclude that the sample means are significantly different. In the other words, the tutoring program developed by the school had significant impact on the SAT score of its students.

We can confirm this by running a formal hypothesis testing with scipy.stats.ttest_rel(). Note that this assumes dependent t-test.

The computed pvalue=0.004 is smaller than the significance level of alpha = 0.05, and therefore we reject the null hypothesis and accept the alternate hypothesis, which is consistent with the conclusion drawn from the confidence interval of difference in mean.

Notes: The above hypothesis testing answers the question of "Did this tutoring program had a significant impact on the SAT scores of students?". However, in cases like this, a more intuitive question is "Did this tutoring program significantly improve the SAT scores of students?" The former uses two-tailed test, and the latter uses one-tailed test, and the procedures for them are a little different.

Checking results with R:

    
        x1 = c(1480, 1280, 890, 340, 1550, 1230, 1010, 1590, 1390, 980)
        x2 = c(1510, 1460, 1320, 700, 1550, 1420, 1340, 1570, 1500, 1300)

        t.test(x1, x2, paired=TRUE)  

        #     Paired t-test

        # data:  x1 and x2
        # t = -3.7753, df = 9, p-value = 0.004381
        # alternative hypothesis: true difference in means is not equal to 0
        # 95 percent confidence interval:
        #  -308.64568  -77.35432
        # sample estimates:
        # mean of the differences 
        #                    -193 
    

4.3. Confidence interval of variance¶

Confidence interval of variance is used to estimate the population variance from sample data and quantify the related uncertainty. C.I. of variance is seldom used by itself, but rather used in conjunction with f-test, which tests equality of variances of different populations. Similar to how the confidence interval of difference in mean forms the foundation of t-test, C.I. of variance forms the foundation of f-test. In the field of statistics and machine learning, the equality of variance is an important assumption when choosing which technique to use. For example, when comparing the means of two samples, student's t-test should not be used when you have a reason to believe that the two samples have different variances. Personally, I found f-test to be useful for the purpose of reading and understanding scientific papers, as many of the papers I have read use f-test to test their hypothesis, or use a variation of f-test for more advanced techniques. It is a pre-requisite knowledge you need to know to understand the more advanced techniques.

I mentioned that different statistics exhibit different distributions above. When a sample data set originates from a normal distribution, its sample means are normally distributed as shown in figure (7). On the other hand, its sample variances are chi-square ($\chi^2$) distributed as shown in figure (10) The curve is asymptotic, and never touches the x-axis. The cumulative probabilty, which is often referred to as "p-value" in hypothesis testing, propagates from the right (p-value=0) to the left (p-value=1). For example, $\chi^2_{.975}=2.70$ is in the lower/left-tail and $\chi^2_{.025} = 19.02$ is in the upper/right-tail. When the samples follow a normal distribution, the $\chi^2$ statistic values can be plugged into eq (10) to compute the confidence interval of variance.

Figure 10: 95% confidence interval of variance.

In confidence interval of variance, the degrees of freedom is:

Recall that the goal of any confidence interval is to estimate the population parameter from a fraction of its samples due to the high cost of obtaining measurement data of the entire data set, as explained in Population vs Samples. You attempt to estimate the population variance $\sigma^2$ within the range of uncertainty with the sample variance $s^2$ obtained from a set of n samples that are "hopefully" representative of the true population.

Confidence interval of variance assumes normality of samples, and is very sensitive to the sample distribution's deviation from normality. In case of non-normal sample distributions, you can either 1) transform the distribution to normal distribution with Box-Cox transformation, or 2) use non-parametric alternatives. For practitioners, I do not recommend 1) unless you really understand what you are doing, as the back transformation process of Box-Cox transformation can be tricky. Furthermore, it doesn't always result in successful transformation of non-normal to normal distribution, as discussed below. I recommend to use 2). If you have non-normal samples and your goal is to compute the C.I. of variance, use bootstrap. If your goal is to check the equality of variances of multiple sample data sets with hypothesis testing, use Levene's test. Both are the non-parametric alternatives that does not require normality of samples.

Pythonic Tip: Computing confidence interval of variance

Unfortunately, there's no Python or R library that computes the confidence interval of variance. The fact that the pre-built function does not exist both in Python and R suggests that the C.I. of variance is seldom used. But as I mentioned before, the reason that I introduce the C.I. of variance is because it forms the foundation of f-test, a statistical hypothesis test that is widely used in scientific papers. The C.I. of variance can be manually computed with eq (10). Don't forget to compute sample variance, instead of population variance by setting ddof=1 as explained above.

The output suggests that the 95% confidence interval of variance is — $\text{C.I.}_{variance}: \,\, 0.072 < \sigma^2 < 0.582$

4.4. Confidence interval of other statistics: Bootstrap¶

(Note: For those people who have web-development experience, this is not CSS Bootstrap.)

I mentioned that different formulas are used to construct confidence intervals of different statistics above. There are three problems with computing the confidence interval of statistics with analytical solutions:

1. Not all statistics have formulas for their confidence intervals

2. Their formulas can be so convoluted, that it may be better to use numerical alternatives

3. You have to memorize their formulas

Bootstrapping is nice because it allows you to avoid these practical concerns. For example, there are no formulas to compute the confidence interval of covariance and median. On the other hand, regression coefficient has its own formula for its confidence interval, but the formulas get really messy in cases of multi-linear or non-linear regression. Wouldn't it be nice if there's a "magic" that saves you from all the math you have to worry about?

Bootstrapping is a statistical method for estimating the sampling distribution of a statistic by sampling with replacement from the original sample, most often with the purpose of estimating confidence intervals of a population parameter like a mean, median, proportion, correlation coefficient or regression coefficient.

Bootstrap can construct confidence intervals of any statistics when combined with Monte Carlo method. The process is visually shown in fig (?). Initially you have 5 samples $[8, 5, 4, 6, 2]$ that you collected from an unknown population. You randomly draw $n=5$ samples from the original sample pool WITH REPLACEMENT, and they become your single bootstrap sample. You repeat this process $r=6$ times to collect multiple bootstrap samples. For each bootstrap sample, you run your functions to compute the statistic of your interest: in this case, np.mean(single_boot). Now you have $r=6$ sample means obtained from $r$ bootstrap samples. You can construct 95% confidence interval of mean with percentile method: np.percentile(mutliple_boot_means, 97.5), np.percentile(mutliple_boot_means, 2.5)

Note that the bootstraped samples will contain duplicate elements a lot, due to random sampling WITH REPLACEMENT. This causes problems with bootstrapping regression models, as explained below.

Figure 11: Bootstrap 95% confidence interval of mean.

Why does bootstrapping work?

(Before you read this section, make sure you understand the difference between Population vs Samples.)

Practitioners wonder WHY bootstrapping works: why is it that resampling the same sample over and over gives good results? If we are resampling from our sample, how is it that we are learning something about the population rather than only about the sample? There seems to be a leap which is somewhat counter-intuitive. The idea comes from the assumption that the sample is a reasonable representation of its underlying population — the population is to the sample as the sample is to the bootstrap samples.

You want to ask question of a population, but you can't because you lack the resources to get measurement data of all possible data points. So you take a fraction of the population, a.k.a the sample, and ask the question of it instead. Now, how confident you should be that the sample answer is close to the population answer depends on how well the sample represents the underlying population. One way you might learn about this is to take samples from different portions of the population again and again. You ask the same question to the multiple samples you collected, and see the variability of the different sample answers to quantify the related uncertainty of your estimation. When this is not possible due to practical limitations, you make some assumptions about the population (ex: population is normally distributed), or use the information in the collected sample to learn about the population.

In bootstrapping, you treat the original sample (size=$n$) you randomly acquired from the population as if it's the population itself of size $n$. From the original sample that you treated to be the population, you randomly draw samples, each of size $n$, with replacement multiple times to simulate direct sampling from the original population. By doing so, you essentially imitate sampling different portions of the original population multiple times. This idea is shown in figure (12).

This is a reasonable thing to do for two reasons. First, the sample in your had is the best you've got, indeed the only informatoin you have about the population. Second, if the original sample is randomly chosen, it will look like the original population they came from. This means that the sample is a good representation of its underlying population. However, if its not a good representation of the population, bootstrap fails. In fact, there's not much you can do in the first place if your sample is biased.

Figure 12: Intuitive idea behind Bootstrapping

Assumptions and Limitations of Bootstrap

Bootstrapping is great because it saves you from the normality assumption of distributions and all the math you have to know to construct confidence intervals. However, just like many other techniques, bootstrap has its own caveats. While bootstrap is distribution-free, it is not assumption-free. The assumptions are listed in this section.

Please note that there is a humongous variety of the bootstrap procedures, each addressing the particular quirk in either the statistic, the sample size, the dependence, or whatever an issue with the bootstrap could be. I am not introducing all of them here as the in-depth technical discussion of bootstrap needs another devoted post, but I still want you to know some of the critical assumptions; I want you to know what you don't know, so that you can google later to learn in-depth.

1. A sample is a good representation of its underlying population

1. Insufficient samples make the bootstrap C.I. to be narrower than the analytical C.I.

There's a myth in the field of statistics that bootstrap is a "cure" for small sample size. NO, it's not. First, if you have too small sample, by a high chance it is not diverse enough to represent all (reasonably) possible aspects of its population. Therefore, it is not a good representation of its population. Second, small sample size makes its bootstrap C.I. to be narrower than the analytical C.I.. This means that bootstrap C.I. reports small uncertainty even when the sample size is small. Not only this is counter-intuitive, but also it is a violation of the mathematic property of C.I. described by eq (1); small sample size $n$ in the denominator of eq (1) should give wider C.I.. But this is not true with bootstrap C.I. as shown in the below simulation result in figure (13).

Three things to note in the figure. First, the upper & lower error bars of bootstrap C.I. of means are asymmetric. This is because bootstrap C.I. is not based on $\pm$ standard error method. This is a very useful property to estimate the central tendency of asymmetric (skewed) populations. Second, both bootstrap and analytical C.I. become narrower with the increasing sample size. This intuitively and mathematically makes sense. Third, bootstap C.I. approximates the analytical C.I. very well with large sample size. This is perhaps the most important advantage of using bootstrap. If you have large sample size, you really don't have to worry anything else (except the indepence of samples), and just stick to bootstrap. All the disadvantages of bootstrap will be overcome by the large sample size.

One might wonder what is "large" enough in practical applications. Unfortunately, the definition of "large" is different for every applications. In the simulation result of figure (13), it seems that $n = 20$ falls in the category of "large" to approximate C.I. of the mean of a normally distributed population with bootstrap. On the other hand, $n=100$ seems to be "large" in case of C.I. of the variances of a normally distributed population. The definition of "large $n$" can vary with different applications (ex: non-normal data, C.I. of regression coefficient or covariance). Carefully investigate your samples to have a good definition of "large".

Figure 13: Comparison of Bootstrap vs Analytical C.I. for different sample sizes

Source Code For Figure (13)

                
                    import matplotlib.pyplot as plt
                    import numpy as np
                    from scipy import stats

                    sample_sizes = [3, 5, 8, 10, 15, 20, 30, 50, 80, 100, 500, 1000]
                    mean = 4            # sample mean
                    std = 3             # sample standard deviation
                    boot_iter = 10000   # bootstrap iterations

                    boot_mean_lo = np.array([])
                    boot_mean_hi = np.array([])
                    analy_mean_lo = np.array([])
                    analy_mean_hi = np.array([])

                    boot_var_lo = np.array([])
                    boot_var_hi = np.array([])
                    analy_var_lo = np.array([])
                    analy_var_hi = np.array([])

                    means = np.array([])
                    variances = np.array([])

                    for size in sample_sizes:

                        np.random.seed(size * 5)

                        arr = np.random.normal(mean, std, size) # randomly draw from a normal distribution

                        # analytical confidence interval of mean
                        means = np.append(means, np.mean(arr))
                        analy_conf_mean = stats.t.interval(1 - 0.05, len(arr) - 1, loc=np.mean(arr), scale=stats.sem(arr))
                        analy_mean_lo = np.append(analy_mean_lo, analy_conf_mean[0])
                        analy_mean_hi = np.append(analy_mean_hi, analy_conf_mean[1])

                        # bootstrap confidence interval of mean
                        boot_means = [np.mean(np.random.choice(arr, len(arr))) for _ in range(boot_iter)]
                        boot_mean_lo = np.append(boot_mean_lo, np.percentile(boot_means, 2.5))
                        boot_mean_hi = np.append(boot_mean_hi, np.percentile(boot_means, 97.5))

                        # analytical confidence interval of variance
                        variances = np.append(variances, np.var(arr, ddof=1))
                        analy_conf_var = (
                            (len(arr) - 1) * np.var(arr, ddof=1) / stats.chi2.ppf(1 - 0.05 / 2, len(arr) - 1),
                            (len(arr) - 1) * np.var(arr, ddof=1) / stats.chi2.ppf(0.05 / 2, len(arr) - 1)
                        )
                        analy_var_lo = np.append(analy_var_lo, analy_conf_var[0])
                        analy_var_hi = np.append(analy_var_hi, analy_conf_var[1])

                        # bootstrap confidence interval of variance
                        boot_vars = [np.var(np.random.choice(arr, len(arr)), ddof=1) for _ in range(boot_iter)]
                        boot_var_lo = np.append(boot_var_lo, np.percentile(boot_vars, 2.5))
                        boot_var_hi = np.append(boot_var_hi, np.percentile(boot_vars, 97.5))


                    # plotting

                    def styling(ax, xticks, xticklables):
                        ax.legend(fontsize=14, loc='lower right', framealpha=1, frameon=True)
                        ax.set_xlabel('Sample sizes', fontsize=16)
                        ax.set_facecolor('#eeeeee')
                        ax.grid(True, linestyle='--', color='#acacac')
                        ax.tick_params(color='grey')
                        ax.set_xticks(xticks)
                        ax.set_xticklabels([str(label) for label in xticklables])
                        _ = [spine.set_edgecolor('grey') for spine in ax.spines.values()]


                    x = np.array([i for i in range(len(sample_sizes))])
                    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

                    axes[0].errorbar(x - 0.15, means, yerr=[abs(boot_mean_lo - means), abs(boot_mean_hi - means)],
                                fmt='o', label='95% bootstrap C.I.', color='k', markersize=8, capsize=5, linewidth=2)
                    axes[0].errorbar(x + 0.15, means, yerr=np.array([abs(analy_mean_hi - means), abs(analy_mean_lo - means)]),
                                fmt='o', label='95% analytical C.I.', color='grey', markersize=8, capsize=5, linewidth=2)

                    styling(axes[0], x, sample_sizes)
                    axes[0].set_ylabel('Sample mean', fontsize=16)
                    axes[0].set_title('Confidence interval of means $\mu$', fontsize=18)
                    axes[0].text(0.75, 0.85, 'aegis4048.github.io', fontsize=15, ha='center', va='center',
                            transform=axes[0].transAxes, color='grey', alpha=0.5);

                    axes[1].errorbar(x - 0.15, means, yerr=[abs(boot_var_lo - variances), abs(boot_var_hi - variances)],
                                fmt='o', label='95% bootstrap C.I.', color='k', markersize=8, capsize=5, linewidth=2)
                    axes[1].errorbar(x + 0.15, means, yerr=np.array([abs(analy_var_hi - variances), abs(analy_var_lo - variances)]),
                                fmt='o', label='95% analytical C.I.', color='grey', markersize=8, capsize=5, linewidth=2)

                    styling(axes[1], x, sample_sizes)
                    axes[1].set_ylabel('Sample variance', fontsize=16)
                    axes[1].set_title('Confidence interval of variances $\sigma^2$', fontsize=18)
                    axes[1].text(0.75, 0.35, 'aegis4048.github.io', fontsize=15, ha='center', va='center',
                            transform=axes[1].transAxes, color='grey', alpha=0.5);

                    fig.tight_layout()
                
            

1. Bootstrap fails to estimate extreme quantiles

1. Samples are independent and identically distributed (i.i.d.)

1. Bootstrap iteration (Monte-Carlo method) should be sufficient to reproduce consistent C.I's.

Because bootstrap relies on "random" resampling, the result of any statistical analysis performed with bootstrap can vary from time to time. The extent of variability depends on the number of bootstrap samples $r$, and $r$ should be large enough to guarantee convergence of bootstrap statistics to a stable value. Note that there's a distinction between the size of the original sample $n$ and the number of bootstrap samples $r$. We can't change $n$, but we can change $r$ because $r$ is equivalent to the number of Monte-Carlo iterations, which can be set by a statistician.

So how do we determine what value of $r$ is "large" enough to guarantee convergence of bootstrap statistics? You can do it by obtaining multiple bootstrap analysis results for increasing number of simulations $r$, and see if the result converges to certain range of values, as shown in figure (15). In the figure, it seems that $r=10,000$ is a good choice. But in practice, you want $r$ to be as large as possible, to an extent where the computational cost is not too huge. I've read a research paper where the authors used $r = 500,000$ to really ensure convergence.

Figure 14: Convergence of Bootstrap C.I.

Source Code For Figure (14)

                
                    import matplotlib.pyplot as plt
                    import numpy as np
                    from scipy import stats
                    import pandas as pd


                    # r, or number of bootstrap samples, or number of Monte-Carlo iterations
                    r_boots = [10, 20, 50, 70, 100, 150, 300, 500, 700, 1000, 2000, 10000, 20000, 50000, 100000]

                    size = 50            # original sample size
                    mean = 50            # sample mean
                    std = 13             # sample standard deviation

                    # randomly draw from a normal distribution
                    arr = np.random.normal(mean, std, size)

                    boot_mean_lo = np.array([])
                    boot_mean_hi = np.array([])
                    results = []

                    for r_boot in r_boots:

                        np.random.seed(r_boot)

                        boot_means = [np.mean(np.random.choice(arr, len(arr))) for _ in range(r_boot)]
                        results.append(boot_means)


                    # plotting
                    fig, ax = plt.subplots(figsize=(8, 4))
                    ax.boxplot(results, sym='', whis=[2.5, 97.5], showfliers=False,
                                         boxprops=dict(linewidth=2.0, color='#4e98c3'),
                                         whiskerprops=dict(linewidth=2.0, color='#4e98c3', linestyle='--'), vert=True,
                                         capprops=dict(linewidth=2.0, color='k'),
                                         medianprops=dict(linewidth=2.0, color='#ad203e'))
                    ax.set_title('Convergence of 95% Bootstrap C.I. with increasing Monte-Carlo iterations', fontsize=15)
                    ax.set_ylabel('Sample mean', fontsize=15)
                    ax.set_xlabel('# of bootstrap samples (Monte-Carlo iterations)', fontsize=15)
                    ax.set_xticklabels([str(r_boot) for r_boot in r_boots], rotation=45)
                    ax.set_ylim(41.7, 52.3)
                    ax.set_facecolor('#eeeeee')
                    ax.grid(True, linestyle='--', color='#acacac')
                    ax.tick_params(color='grey')
                    _ = [spine.set_edgecolor('grey') for spine in ax.spines.values()]
                    ax.text(0.21, 0.1, 'aegis4048.github.io', fontsize=15, ha='center', va='center',
                            transform=ax.transAxes, color='grey', alpha=0.5);
                
            

1. Coverage of naive bootstrap is relatively weak compared to more robust bootstrap methods

1. Bootstrapping continuous data is a bit tricky

1. Bootstrap is not robust in heavy-tailed distributions

Pythonic Tip: Bootstrapping in Python

Basic bootstrap

Bootstrap by itself means resampling from a sample. In its simplest form, it can be implemented in just one-line code with np.random.choice. For demonstration, assume that the original sample of size n=500 was randomly drawn from a normal distribution.

Basic bootstrap with Monte-Carlo method

Bootstrap is often combined with Monte-Carlo method to quantify uncertainty in statistics (ex: mean, median, variance, etc...). It just means that you generate multiple instances of the single_boot object above.

With the above simulation, we stored r=1000 instances of bootstrap samples in the object monte_boot. It is a two dimensional array of size $r$ x $n$ (1000 x 500). I will show only the first 10 bootstrap samples with the Pandas DataFrame, since it will be too lengthy if I output all 1000 rows.

Let's visualize the results in Matplotlib to understand Monte-Carlo simulations applied with bootstrap more intuitively. Recall that the goal of Monte-Carlo method is to simulate hundreds or thousands of "possible outcomes". Each line in the right plot below represents a single possible outcome.

Figure 15: 1000 Bootstrap simulations

Basic bootstrap with Monte-Carlo method + constructing confidence intervals

We then analyze the simulation result to do whatever statistical estimation we want to do. Since this article is about confidence intervals, I will show how to construct confidence intervals of various statistics with bootstrap percentile method. It's actually very simple to implement. Just wrap your single bootstrap sample with a function that calculates the statistic of your interest. The result is an array of statistics of 1000 sample sets. In the other words, thousand data points of the statistics.

If the statistic of your interest does not have a library function, you can define your own function, and wrap the bootstrap sample with it.

Let's visualize the 95% confidence interval of various statistics obtained from Monte-Carlo bootstrap. Note that statistics like median, skew, kurtosis do not have analytical solutions to construct C.I., and can only be constructed from numerical methods like Monte-Carlo bootstrap (which makes bootstrap very powerful).

Figure 16: Uncertainty models obtained from Monte-Carlo bootstrap

5. Confidence interval of non-normal distribution¶

One of the biggest assumptions in the field of statistics is the assumption of normality. It is important because of the two key characteristics of normal distribution: symmetry (low skew) and light/short-tails (low kurtosis). Symmetry is important because many statistical techniques use $\pm$ standard error methods that assume equal properties from both sides of a distribution. On the other hand, heavy-tailedness is detrimental in inferencing population parameters because heavy tailed distributions (ex: log-normal distribution, cauchy distribution) have non-negligible number of extreme data points that may introduce bias in your samples (understanding the impact of heavy-tailedness is very important in inferencial statistics. Make sure you understand the issues of it by reading below).

As I mentioned in the Key Takeaway 4 & 8 above, constructing confidence intervals are very different when data is not normally distributed. This section explains the reasons why they are different, and their non-parametric alternatives.

5.1 Problems of non-normal distributions and central tendency¶

When reporting a summary of a sample data, statisticians tend to report the mean (average) of a sample. But have you wondered why they bother specifically about the means? Often times the ultimate goal is NOT to compute a mean of a distribution, but to compute a measure of central tendency of a distribution. A measure of central tendency is a single value that attempts to describe a set of data by identifying the central position within that set of data.

The mean is the measure of central tendency that you are the most familiar with, but there are others, such as the median and the mode. However, mean is not a good measure of central tendency when there is a sign of deviation from normality, which can be characterized by skewness (asymmetry) and kutrosis (heavy-tails).

Problems with Skewness (asymmetry)

Figure (17) below shows the effectiveness of mean, median, and mode as a measure of central tendency for different skewness. For normal distribution (a), most samples are found at the middle, and they are symmetrically distributed around the mean. This makes the mean as a good measure of central tendency. However, in asymmetrical distributions like (b) and (c), the median (or arguably the mode) is a better choice of central tendency, as it is closer to the central location of the distribution than the mean does. Naively reporting the mean as a summary of data results in misinterpretation of the central location of population.

Note that some parametric techniques that assume normality of data are robust to mild skewness. There are many parametric alternatives that account for skewness, such as skewness-adjusted t-test and Box-Cox transformation. More discussion of robustness of parametric techinques relative to skewness is done below.

Figure 17: Central tendency of distributions

Problems with kurtosis (heavy-tails)

(Note that heavy-, long-, fat-tails are all synonyms. Some literature distinguishes among their meanings, but in general, they are all equivalent terms.)

Unlike skewness, which can be relatively easily adjusted, heavy-tails can be very painful in inferencial statistics. When a distribution has heavy-tails, it means that the sample has extreme data points that can introduce bias to statistical inferences (note that the term bias is defined as the difference between the statistical estimation vs. the actual population parameter). This is because the occurrence of extreme data points are so low that they are unlikely to be included in the sample you collected, and yet there are substantial number of them lurking in the uncollected portion of a population that impacts population parameters due to their extremity.

Figure (18) below shows the comparison of heavy- vs light-tailed distributions in a form of probability density function (PDF). Let's talk about the left plot first. Observe that in the heavy-tailed distribution (black line), the chance of occurrence for extreme values does not decay to zero as fast as that of the light-tailed distribution (grey line). This means that you have non-negligible chance of observing extreme values, and they may not have upper/lower bounds. In finance, such property is translated into high risk & unpredictability, which is modeled by Cauchy distributions.

The right plot in figure (18) compares normal vs Cauchy distributions. Note that the tails of Cauchy distribution is higher than that of the normal distributions. The heavy tails in Cauchy attempt to model financial/investment risks by representing randomness in which extremes are encountered. When data arise from an underlying heavy-tailed distribution, shoehorning in the "normal distribution" model of risk would severely understate the true degree of risk. The narrow peak of the Cauchy distribution signifies that less samples fall within the range of Six-Sigma ($\pm 3\sigma$), increasing unpredictability as less samples have the typical values. Practitioners often neglect the distinction between light- vs heavy-tail because they visually don't look very different when plotted in PDF. This is a big misconception.

Figure 18: Heavy-tailed distributions have more extreme data points that can be detrimental to statistical inferences

Treatment method for heavy-tailed distribution depends on the situation. The first step would be to map the distribution to a normal distribution using transformation techniques discussed below. However, the heavy-tails may still be present even after the presentation. If you are running a regression model, you may simply remove the tails and see if there's any performance improvement. In case of inferencial statistics, naive removal maybe detrimental, as it completely ignores the extremity in your population. Sometimes the heavy-tails may suggest amalgamation of two distinct populations, in which case you may need to separate them apart before running your model.

This can be seen easily by running a quick test. I generated two samples from two distinct normal distribution as pop_1 and pop_2, and plotted their combined distributions. Even though the samples were both generated from normal distributions, it exhibit heavy tails due to their difference in central tendencies.

Generating two distinct normal samples:

    
        pop_1 = np.random.normal(10, 3, 100)
        pop_2 = np.random.normal(15, 6, 100)
    

Figure 19: Heavy-tails due to two distinct populations

            
                %matplotlib inline

                import numpy as np
                import matplotlib.pyplot as plt
                import matplotlib.gridspec as gridspec

                np.random.seed(4)

                # generate samples from two distinct populations
                pop_1 = np.random.normal(10, 3, 100)
                pop_2 = np.random.normal(15, 6, 100)

                # combine
                data = np.append(pop_1, pop_2)

                fig = plt.figure(figsize=(10, 4))
                plt.suptitle('Heavy-tails due to amalgamation of two distinct populations', fontsize=15)

                gs = gridspec.GridSpec(2, 2)
                axes = [plt.subplot(gs[0, 0]), plt.subplot(gs[1, 0]), plt.subplot(gs[:2, 1])]

                axes[0].hist(pop_1, histtype='stepfilled', edgecolor='k', alpha=0.4, color='grey', density=True)
                axes[1].hist(pop_2, histtype='stepfilled', edgecolor='k', alpha=0.4, color='grey', density=True)
                axes[2].hist(data, histtype='stepfilled', edgecolor='k', alpha=0.4, color='grey', density=True)

                for ax in axes:
                    ax.text(0.73, 0.85, 'aegis4048.github.io', fontsize=11, ha='center', va='center',
                            transform=ax.transAxes, color='grey', alpha=0.5)
                    ax.grid(True, linestyle='--', color='#acacac')
                    ax.set_xlim(0, 35)

                plt.subplots_adjust(wspace=0.5)
                axes[2].text(-0.415, 0.55, 'Combine', transform=axes[2].transAxes, fontsize=13)
                axes[2].text(-0.4, 0.43, r'$\Longrightarrow$', transform=axes[2].transAxes, fontsize=30)
            
        

Problems with multi-modal distributions

In case of multi-modal distributions, reporting the mean, median, or mode can be uninformative at all, as none of them describes the central location of populations. While it is possible to compute the confidence interval of statistics using non-parameteric methods like bootstrap, it is useless. I hope that figure (20) below is descriptive enough to understand the problems of multi-modal distributions. You will need to devise your own method to describe the central location of data, such as separating the multi-modes into distinct unimodal distributions.

Figure 20: Absence of central tendency in multi-modal distributions

            
                import numpy as np
                import matplotlib.pyplot as plt

                np.random.seed(4)

                n = 1000
                np.random.seed(1)
                x = np.concatenate((np.random.normal(0, 1, int(0.3 * n)),
                                    np.random.normal(14, 1, int(0.7 * n)),
                                    np.random.normal(8, 2, int(0.7 * n))))[:, np.newaxis]

                plt.rcParams.update(plt.rcParamsDefault)

                fig, ax = plt.subplots(figsize=(8, 4))

                adjust = 0.12
                n, bins, patches = ax.hist(x, bins=30, histtype='stepfilled', edgecolor='k', alpha=0.4, color='grey', density=True)
                ax.axvline(np.mean(x) + adjust, color='k', label='Mean', linestyle='-', alpha=0.9)
                ax.axvline(np.median(x) + adjust, color='k', label='Median', linestyle='--', alpha=0.9)
                ax.axvline(bins[np.where(n == max(n))[0][0]] + adjust, color='k',  label='Mode', linestyle=':', alpha=0.9)

                ax.text(0.2, 0.9, 'aegis4048.github.io', fontsize=15, ha='center', va='center',
                        transform=ax.transAxes, color='grey', alpha=0.5)
                ax.grid(True, linestyle='--', color='#acacac')
                ax.set_xlabel('X', fontsize=16)
                ax.set_ylabel('Probability', fontsize=16)
                ax.legend(loc='upper right', framealpha=1, frameon=True)

                fig.tight_layout()
            
        

5.2 Robustness of confidence intervals to non-normality¶

Confidence interval of mean is robust to mild deviation from normality. This means that hypothesis tests that rely on the means (ex: t-test) are also robust to mild deviation from non-normality. In fact, there are a few variations of t-test that deals with mild non-normality (ex: skewness-adjusted t-statistic).

Confidence interval of variance is VERY sensitive to even small deviation from non-normality. Its low popularity in real-life applications can be attributed to its high sensitivity. In real-life, people tend to not care much about the C.I. of variance, but techniques that check the equality of variances of multiple samples with hypothesis testing. When sample is normal, you can use f-test or Barlett's test to check equality of variances. When the sample is not normal, use non-parametric alternative like Levene's test.

Many techniques that assume normality of data are still resistant to mild departure from its assumptions. Then the natural question is, how do I know the severity of deviation from normality? I found from here that the skewness and kurtosis of $-2 < 0 < +2$ is an acceptable deviation from normality, but I've seen some people arguing $-1 < 0 < +1$ is the limit. This all depends on how much you are willing to tolerate the deviation. Another way to assess departure from normality is visualizing the distribution in Q-Q plots, as shown in figure (21). You can also use this online interactive QQ-plot generator to study the effect of asymmetry and heavy-tailedness. Note that I introduced the code snippet for visualizing Q-Q plots in Python below.

Figure 21: Assessing deviation from normality with Q-Q plots

Pythonic Tip: Q-Q plots with SciPy

Q-Q plots are useful when assessing deviation from a pre-defined distribution. It is implemented in SciPy pakcage as scipy.stats.probplot. Note that the pre-defined distribution can be anything of your choice by changing the argument dist; it can be lognormal, Weibull, exponential, Cauchy, or anything. The red line in the Q-Q plot represents the theoretical distribution of your data in a form of Cumulative Density Function (CDF).

5.3 Transform to normal distribution: Box-Cox¶

When you want to construct a confidence interval of a sample, but the sample is non-normal, you have two options: transform the sample into a normal distribution, or use non-parametric alternatives. This section discusses the transformation part.

There are three techniques used to map data to a normal distribution: Box-Cox, Yeo-Johnson, and normal quantile transform (check the SciPy documentation to see the comparison plot of each tecnnique). I will discuss Box-Cox transformation here.

Box-Cox transformation is a type of power transformation to convert non-normal data to normal data by raising the distribution to a power of lambda ($\lambda$). The algorithm automatically solves for $\lambda$ that best transforms a distribution into a normal distribution.

Box-Cox transformation is a statistical technique known to have remedial effects on highly skewed data. If $\lambda$ is determined to be 2, then the distribution will be raised to a power of 2 — $Y^2$. The exception to this rule is when the $\lambda$ is 0 - log will be taken to the distribution — log($Y$). Note that using a $\lambda = 1$ does not do anything to the distribution. If the Box-Cox algorithm spits out $\lambda = 1$, it probably means that your data is Gaussian-like or Gaussian enough to an extent that there is no need for transformation. The figure below shows the effect of raising the exponential distribution (left) with the power of $\lambda = .316$ to acquire the transformed normal distribution (right).

Figure 22: Box-cox transformation

            
                from scipy import stats
                import matplotlib.plyplot as plt

                # generate sample data
                x = stats.expon.rvs(size=100, random_state=1)

                # Box-Cox transform
                xt, lmbda = stats.boxcox(x)

                # plot
                fig, axes = plt.subplots(2, 2, figsize=(10, 5))

                axes[0, 0].hist(x, histtype='stepfilled', edgecolor='k', alpha=0.4, color='grey', density=True)
                axes[0, 1].hist(xt, histtype='stepfilled', edgecolor='k', alpha=0.4, color='grey', density=True)
                stats.probplot(x, dist=stats.norm, plot=axes[1, 0])
                stats.probplot(xt, dist=stats.norm, plot=axes[1, 1])

                axes[0, 0].set_title('Original data')
                axes[0, 1].set_title(r'Transformed data,  $\lambda$=%.3f' % lmbda)
                axes[1, 0].set_title('Normal Q-Q plot')
                axes[1, 1].set_title('Normal Q-Q plot')
                axes[1, 1].set_ylabel('')

                for ax in axes.flatten():
                    ax.grid(True, linestyle='--', color='#acacac')
                    ax.text(0.05, 0.85, 'aegis4048.github.io', fontsize=12, transform=ax.transAxes, color='grey', alpha=0.5)

                for ax in axes.flatten()[2:]:
                    ax.get_lines()[0].set_marker('o')
                    ax.get_lines()[0].set_color('grey')
                    ax.get_lines()[0].set_markeredgecolor('black')
                    ax.get_lines()[0].set_markersize(5.0)
                    ax.get_lines()[0].set_alpha(0.5)
                    ax.get_lines()[1].set_linewidth(1.0)
                    ax.get_lines()[1].set_color('k')

                plt.subplots_adjust(wspace=13)
                fig.text(0.47, 0.33, 'Box-Cox', size=14, ha='center', va='top')
                fig.text(0.465, 0.3, r'$\Longleftrightarrow$', size=30, ha='center', va='top')
                fig.text(0.47, 0.8, 'Box-Cox', size=14, ha='center', va='top')
                fig.text(0.465, 0.77, r'$\Longleftrightarrow$', size=30, ha='center', va='top')
                fig.tight_layout()
            
        

Not all transformations work perfectly

The example illustrated above is the most ideal case of Box-Cox transform. Not all transforms work perfectly, especially in case of heavy tails. Visualize your transformed histrogram, and the Q-Q plot to evaluate the performance of the transform. If the transform did not mitigate the skewness (this is rare), you will have to look for non-parametric alternatives, like bootstrap explained above. If the transform did not mitigate the heavy-tails, your solution depends on the situation, as discussed above.

Figure 23: Box-cox transformation

            
                from scipy import stats
                import matplotlib.pyplot as plt

                # generate sample data
                x = stats.skewnorm.rvs(100, size=1000) + 1

                # Box-Cox transform
                xt, lmbda = stats.boxcox(x)

                # plot
                fig, axes = plt.subplots(2, 2, figsize=(10, 5))

                axes[0, 0].hist(x, histtype='stepfilled', edgecolor='k', alpha=0.4, color='grey', density=True)
                axes[0, 1].hist(xt, histtype='stepfilled', edgecolor='k', alpha=0.4, color='grey', density=True)
                stats.probplot(x, dist=stats.norm, plot=axes[1, 0])
                stats.probplot(xt, dist=stats.norm, plot=axes[1, 1])

                axes[0, 0].set_title('Original data')
                axes[0, 1].set_title(r'Transformed data,  $\lambda$=%.3f' % lmbda)
                axes[1, 0].set_title('Normal Q-Q plot')
                axes[1, 1].set_title('Normal Q-Q plot')
                axes[1, 1].set_ylabel('')

                for ax in axes.flatten():
                    ax.grid(True, linestyle='--', color='#acacac')
                    ax.text(0.05, 0.85, 'aegis4048.github.io', fontsize=12, transform=ax.transAxes, color='grey', alpha=0.5)

                for ax in axes.flatten()[2:]:
                    ax.get_lines()[0].set_marker('o')
                    ax.get_lines()[0].set_color('grey')
                    ax.get_lines()[0].set_markeredgecolor('black')
                    ax.get_lines()[0].set_markersize(5.0)
                    ax.get_lines()[0].set_alpha(0.5)
                    ax.get_lines()[1].set_linewidth(1.0)
                    ax.get_lines()[1].set_color('k')

                plt.subplots_adjust(wspace=13)
                fig.text(0.47, 0.33, 'Box-Cox', size=14, ha='center', va='top')
                fig.text(0.465, 0.3, r'$\Longleftrightarrow$', size=30, ha='center', va='top')
                fig.text(0.47, 0.8, 'Box-Cox', size=14, ha='center', va='top')
                fig.text(0.465, 0.77, r'$\Longleftrightarrow$', size=30, ha='center', va='top')
                fig.tight_layout()
            
        

All data to be positive and greater than 0 (Y > 0)

Box-Cox transformation does not work if data is smaller than 0. This can easily be fixed by adding a constant ($C$) that will make all your data greater than zero. The transformation equation is then:

Note that you do not have to worry about this if you are using Yeo-Johnson transform instead of Box-Cox transform.

Back-transformation

Since box-cox transform raises a sample to the power of $\lambda$, the scale of your sample changes. While removing skewness in a sample is desirable, the change in scale is not. Once you construct the confidence interval of a statistic of your interest, you have to take the inverse of the $\lambda$ power to your estimated statistics to put it back to its original scale. In the same way, if you've added a constant $C$ to satisfy the assumption above, if have to subtract your statistical estimation by $C$ at the end.

Pythonic Tip: Box-Cox transform with SciPy and Scikit-Learn

There are two ways to do Box-Cox transform, and its back transform: with scipy.special.boxcox and sklearn.preprocessing.PowerTransformer. Scipy's implementation is preferrable for one-dimensional data, and Scikit-Learn's for multi-dimensional data. Note that you also need to backtransform your data, or your calculated statistics to its original scale. Here, I will calculate the centeral tendency of non-normal distribution, which is described by a median.

SciPy Implementation

First, you need to transform the data with scipy.special.boxcox. Make sure to return & store lmbda, which raises the original distribution to the power of $\lambda$. You need it later to back-transform the calculated statistic into its original scale.

Let's first generate non-normal sample data we are going to use, and take a look at it. Assess the normality of your sample with Q-Q plots explained above.

The distribution seems to be non-normal. The mean is 8.30, and the median is 5.26. Let's transform this data into a normal distribution with Box-Cox transform, and see how it changes.

The transformation seemed to work well. Note the change in scale due to transform. The mean changed from 8.30 to 1.63, and the median changed from 5.26 to 1.64. This shows the need to apply back-transformation to the calculated statistic.

We proceed to calculate the confidence interval of your statistic. In this case, we compute the 95% C.I. of mean, with the code snippet described above. If the transformation did not work well, use a non-parametric alternative explained below.

Note that the above C.I. of mean is in the transformed scale. We need to apply back-transform to put it back to its original scale. We use scipy.special.inv_boxcox.

Taking the inverse of the transformed C.I. of mean returned a range of 4.93 ~ 5.54. Recall that in the original sample, Mean = 8.30, and median = 5.26. The result is consistent with the statement in the WARNING above: "when the transformed variable is symmetric, taking an inverse of the transformed mean yields the median of the original variable."

Additionally, we can also try a non-parametric alternative, bootstrap, to obtain the C.I. of median of the non-normal distribution. Recall that the variable x generated above is a non-normal sample. One can see that the C.I. of median obtained from Box-Cox and bootstrap agree with each other.

Scikit-Learn implementation

Box-cox transform can also be implemented with Scikiy Learn: sklearn.preprocessing.PowerTransformer. Scikit-Learn's implementation of box-cox transformation is specialized for multi-dimensional data, and therefore it accepts 2-D array as its input. It's ability to handle multi-dimensional data is quite useful in machine learning.

This function is available in versions scikit-learn v0.20 or above. Most versions of Anaconda support scikit-learn v0.19.1. If you are using Anaconda Python, you will need to update Scikit Learn with the following command in the console (note that if you are using Jupyter Notebook, you must restart your Jupyter after upgrading Scikit-Learn):

For consistency, we will keep using the non-normal sample x I generated above. Before we start, recall that x is an 1-D array. Since Scikiy-Learn's Box-Cox expects 2-D array as input, you will have to change the data type into a numpy matrix with np.asmatrix:

(Note: You need to have some basic understanding of Objected Oriented Programming (OOP) to understand this explanation)

Let's apply the transform. First, instantiate the transformer object pt. Second, fit the $\lambda$ parameter using fit(), which will store the $\lambda$ parameter as a class attribute inside the pt object. The $\lambda$ parameter stored will be used later to apply inverse transform. Last, use transform() to apply Box-Cox transformation.

Note that appling pt.fit() and pt.transform() can also be done with one-line with pt.fit_transform().

From this point, the procedures for obtaining the C.I. of median of a non-normal distribution are the same as in SciPy's implementation, except that we use pt.inverse_transform (Scikit-Learn) instead of inv_boxcox (SciPy), and convert 1-D array to 2-D datatype array. Note that the output is the exact same is the one given by the SciPy implementation above.

5.4 Non-parametric alternative: Bootstrap¶

(Note: This section is a continuation of Section 4.4. Confidence interval of other statistics: Bootstrap explained above. For bootstrap code, navigate to here.)

Bootstrap is the most popular choice when making a point estimation of a population parameter of an unknown distribution. It is distribution-free. This means that it does not assume anything about the shape of the distribution (Ex: normal, lognormal, weibull, exponential). For example, computing the confidence interval of mean with eq (1) requires symmetry of a distribution, as it relies on $\pm$ standard error method. Any parametric methods (Ex: t-test, f-test) that assume normality of data has less coverage than their non-parametric alternatives when their assumptions are violated.

Let's take a deeper look at the robustness of parametric vs non-parametric methods for different distributions by running simulations. We compute the confidence interval of mean and variance. CI of mean is known to be resistent to deviation from normality, while CI of variance is very sensitive. The procedure is as follows:

1. Randomly generate POPUPLATION_SIZE observations from a specified distribution. Assume this to be a population. We generate 4 populations: normal, exponential, skewed, and lognormal. Each population goes through step 2-5.

2. Randomly draw samples NUM_SAMPLES times from a population. Each sample has a size = SAMPLE_SIZE. Each sample goes through step 3-4.

3. Compute the parametric & non-parametric 95% confidence interval of variance and mean for a sample. For parametric CI of mean and variance, use eq (1) and eq (10). For non-parametric CI, use bootstrap.

4. Check if a population mean or variance falls within the computed confidence interval.

5. Count the number of times a population mean and variance falls within the computed confidence interval.

The above procedure describes a situation in which we attempt to inference the unknown population parameters from a small portion of the population, a.k.a. the samples (check Section 3. Population vs Samples above). We want to know if the computed confidence interval of samples contain the true population parameter. Figure (24) summarizes the simulation result:

Highlights of the simulation result

There are lots of things we can learn from the simulation result, shown in figure (24). I want to mainly emphasize on the distinction between the parametric vs. non-parametric methods here. Note that bootstrap is one kind of non-parametric methods

1. Some statistics are sensitive to deviation from normality

1. Parametric methods acquire 95% nominal coverage rate ONLY under normality

1. Making an inference about heavy-tailed distribution is difficult.

1. Bootstrap outperforms parametric method under non-normality

1. Parametric method outperforms non-parametric method under normality

1. Non-parametric methods yield "approximate" solutions with smaller risk