This tutorial explains how to build confidence regions (the 2D version of a confidence interval) using as little statistical theory as possible. I also avoid the traditional terminology and notation such as *α*, Z_{1-α}, critical value, confidence level, significance level and so on. These can be confusing to beginners and professionals alike.

Instead, I use simulations and two keywords only: confidence region, and confidence level. The purpose is to explain the concept using a framework that will appeal to machine learning professionals, software engineers and non-statisticians. My hope is that you will gain a deep understanding of the technique, without headaches. I also introduce an alternative type of confidence region, called *dual confidence region*. It is asymptotically equivalent to the standard definition. In my opinion, it is more intuitive.

## Example

This example comes from a real-life application. In this section I provide the minimum amount of material necessary to illustrate the methodology. The full problem is described in the last section, for the curious reader. In its simplest form, we are dealing with independent bivariate Bernoulli trials. The data set has *n* observations. Each observation consists of two measurements (*u _{k}*,

*v*), for

_{k}*k*=1, …,

*n*. Here

*u*= 1 if some interval

_{k}*B*contains zero point (otherwise

_{k}*u*= 0). Likewise,

_{k}*v*= 1 if the same interval contains one point (otherwise

_{k}*v*= 0).

_{k}The interval *B _{k}* can contain more than one point, but of course it can not simultaneously contain one and two points. The probability that

*B*contains zero point is

_{k}*p*; the probability that it contains one point is

*q*, with 0<

*p*+

*q*<1. The goal is to estimate

*p*and

*q*. The estimators (proportions computed on the observations) are denoted as

*p*and

_{0}*q*.

_{0}Since we are dealing with Bernoulli variables, the standard deviations are *σ _{p}* = [

*p*(1-

*p*)]

^{1/2}and

*σ*= [

_{q}*q*(1-

*q*)]

^{1/2}. Also the correlation between the two components of the observation vector is

*ρ*

_{p,q}= -pq /

*σ*. Indeed the probability to observe (0, 0) is 1-

_{p}σ_{q}*p*–

*q*, the probability to observe (1, 0) is

*p*, the probability to observe (0, 1) is

*q*, and the probability to observe (1, 1) is zero.

## Shape of the Confidence Region

A confidence region of level *γ* is a domain of minimum area that contains a proportion *γ* of the potential values of your estimator (*p*_{0}, *q*_{0}), based on your *n* observations. When *n* is large, (*p*_{0}, *q*_{0}) approximately has a bivariate normal distribution (also called Gaussian), thanks to the central limit theorem. The covariance matrix of this normal distribution is specified by *σ _{p}*,

*σ*and

_{q}*ρ*

_{p,q}measured at

*p*=

*p*

_{0}and

*q*=

*q*

_{0}. For a fixed

*γ*, the optimum shape — the one with minimum area — necessarily has a boundary that is a contour level of the distribution in question. In our case, that distribution is bivariate Gaussian, and thus contour levels are ellipses.

Let us define

H_n(x,y,p,q)=\frac{2n}{1-\rho_{p,q}^2}\cdot \Big[\Big( \frac{x-p}{\sigma_p}\Big)^2 -2\rho_{p,q}\Big(\frac{x-p}{\sigma_p}\Big)\Big(\frac{y-q}{\sigma_q}\Big) + \Big(\frac{y-q}{\sigma_q}\Big)^2\Big].

This is the general elliptic form of the contour line. Essentially, it does not depend on *n*, *p*, *q* when *n* is large. The standard confidence region is then the set of all (*x*, *y*) satisfying *H _{n}*(

*x*,

*y*,

*p*

_{0},

*q*

_{0}) ≤

*G*. Here you choose

_{γ}*G*to guarantee that the confidence level is

_{γ}*γ*. Replace ≤ by = to get the boundary of that region.

In this case *G _{γ}* is a quantile of the Hotelling distribution. In the simulation section, I show how to compute

*G*. The simulations apply to any setting, whether

_{γ}*G*is a Hotelling, Fisher or any quantile. Or whether the limit distribution of your estimator (

_{γ}*p*

_{0},

*q*

_{0}) is Gaussian or not, as

*n*— the sample size — increases. These simulations provide a generic framework to compute confidence regions.

## Dual Confidence Region

The dual confidence region is simply obtained by swapping the roles of (*x*, *y*) and (*p*, *q*) in *H _{n}*(

*x*,

*y*,

*p*,

*q*). It is thus defined as the set of (

*x*,

*y*) satisfying

*H*(

_{n}*p*,

*q*,

*x*,

*y*) ≤

*H*. Again, you choose

_{γ}*H*to guarantee that the confidence level is

_{γ}*γ*. Also, (

*p*,

*q*) is replaced by (

*p*

_{0},

*q*

_{0}). This is no longer the equation of an ellipse. In practice, both confidence regions are very similar. Also,

*H*is almost identical to

_{γ}*G*. The interpretation is as follows. A point (

_{γ}*x*,

*y*) is in the dual confidence region of (

*p*

_{0},

*q*

_{0}) if and only if (

*p*

_{0},

*q*

_{0}) is in the standard confidence region of (

*x*,

*y*). We use the same

*n*and confidence level

*γ*for both regions. You can use the same principle to define dual confidence intervals.

Figure 1 shows an example based on simulations.

## Simulations

The simulations consist of generating *N* data sets, each with *n* observations. Use the joint Bernoulli model described in the first section, for the simulations. The purpose is to create data sets that have the same statistical behavior as your observations. In particular, use *p*_{0} and *q*_{0} in the bivariate Bernoulli model, for the simulations.

For each simulated data set, compute the proportions, standard deviations and correlations. They are denoted as *x* , *y*, *σ _{x}*,

*σ*and

_{y}*ρ*

_{x,y}(one set of values per data set). Use the standard formulas from this article: for instance,

*σ*= [

_{x}*x*(1-

*x*)]

^{1/2}. Also compute

*G*(

*x*,

*y*) =

*H*(

_{n}*x*,

*y*,

*p*

_{0},

*q*

_{0}) and

*H*(

*x*,

*y*) =

*H*(

_{n}*p*

_{0},

*q*

_{0},

*x*,

*y*) for each data set. Put the results in a table with

*N*rows and 7 columns. Proceed as follows.

**Standard confidence regions**: sort the table by*G*(*x*,*y*).**Dual confidence region**: sort the table by*H*(*x*,*y*).

The first *γN* rows in your sorted table determines your confidence region of level *γ*. All the (*x*, *y*) in those rows belong to your confidence region. In the first *γN* rows, the last value of *H*(*x*, *y*) — if sorted by *H*(*x*, *y*) — is *H _{γ}*. Likewise, the last value of G(

*x*,

*y*) — if sorted by G(

*x*,

*y*) — is

*G*. See example in Figure 1, with

_{γ}*N*= 10,000 and

*n*= 20,000. As

*N*increases, your simulations yield regions closer and closer to the theoretical ones. The spreadsheet with these simulations is available on my GitHub repository, here.

## The Original Problem

The original problem consisted of estimating the two parameters of a perturbed lattice point process. These stochastic processes have applications in sensor and cell network optimization. Rather than a direct estimation, I used proxy statistics *p*, *q* for the estimator. This method, called minimum contrast estimation, requires a one-to-one-mapping between the original parameter space, and the proxy space.

The point count statistic discussed earlier measures the number of points of this process that are in a specific interval *B _{k}*. I used

*n*non-overlapping intervals

*B*

_{1}, …,

*B*, each one yielding one observation vector. The observation vectors are almost identically and independently distributed across the intervals. However, the first and second components of the vectors are negatively correlated. This explains the choice of the bivariate Bernoulli distribution for the model. The topic is discussed in details in my upcoming book, here.

_{n}#### About the Author

Vincent Granville is a machine learning scientist, author and publisher. He was the co-founder of Data Science Central (acquired by TechTarget) and most recently, founder of Machine Learning Recipes.