Why does this formula produce $p_{2}$ probabilities for Mahalanobis distances?

Question

At the IBM website it is written that

The p1 probabilities are standard probabilities of an observation from a multivariate normal distribution being that far or further from the centroid, which are based on a cumulative central chi-squared distribution with degrees of freedom equal to the number of observed variables.

and

The p2 probabilities are probabilities of the ordered values of N distances being as far or further away from the centroid.

The AMOS development website adds that

Small numbers in the p1 column are to be expected. Small numbers in the p2 column, on the other hand, indicate observations that are improbably far from the centroid under the hypothesis of normality.

The formula for $p_{2}$, along with a worked example is given at another part of the same website. The page discusses “Calculating p2 for the case with the k-th largest d2”, stating that

In general, for the case that is k-th furthest from the centroid (meaning that there are $k-1$ cases further from the centroid), $p_{2}$ is calculated by first evaluating $p_{1}$ for that case and then calculating
$p_{2} = 1 \\ - _{N}C_{N-0}(1-p_{1})^{N}(p_{1})^{0} \\ - _{N}C_{N-1}(1-p_{1})^{N-1}(p_{1})^{1} \\ - _{N}C_{N-2}(1-p_{1})^{N-2}(p_{1})^{2} \\ ...\\ - _{N}C_{N-k+1}(1-p_{1})^{N-k+1}(p_{1})^{k-1} \\ $

where N is the number of cases.

I understand how to convert d-squared values into $p_{1}$ values. With some R code (below) I can reproduce the $p_{2}$ values given in the worked example at the aforementioned webpage. I also understand the point of the $p_{2}$ values is to provide the probability the case from our dataset with this rank (highest d-squared, second-highest d-squared, etc) would be that far or further from the centroid, assuming multivariate normality. So if the highest d-squared values come along with very low $p_{2}$ values, that suggests a violation of multivariate normality.

But I don’t really understand why the formula should give these probabilities. What is the intuition behind it?

amos_p1 <- c(.0046132,.0085718,.0390278,.0437704,.0475222)

N = 73;

amos_p2 <- numeric(length(amos_p1))


for (i in 1:length(amos_p1))
  { 
  k <- i;
  p1_value <- amos_p1[i];

  start_value <- 1;

  while (k >= 1)
    {
      start_value = start_value - choose(N,N-k+1) * (1-p1_value)^(N-k+1) * (p1_value)^(k-1)
      k <- k-1;
    }

  amos_p2[i] <- start_value;

  } 

print(amos_p2)

Answer 1

These values are order statistics, so they follow standard distributional results for those types of statistics. In the case where you have random vectors $\mathbf{X}_1, ..., \mathbf{X}_n \sim \text{IID N}(\boldsymbol{\mu}, \mathbf{\Sigma})$ from the multivariate normal distribution with dimension $m$, you get the squared Mahalanobis distances:

$$D_i \equiv D^2(\mathbf{X}_i) = (\mathbf{X}_i - \boldsymbol{\mu})^\text{T} \mathbf{\Sigma}^{-1} (\mathbf{X}_i - \boldsymbol{\mu}).$$

It can be shown that $D_i \sim \text{ChiSq}(m)$, so if we let $d_i = D^2(\mathbf{x}_i)$ denote the corresponding observed value, then the "$p_1$" probabilities for these observed distances are:

$$\begin{aligned} p_{1}(d_i) &\equiv \mathbb{P}(D_i \geqslant d_i) \\[6pt] &= \int \limits_{d_i}^\infty \text{ChiSq}(r|m) \ dr \\[6pt] &= 1 - F_\text{ChiSq}(d_i|m). \\[6pt] \end{aligned}$$

Now, letting $D_{(1)},...,D_{(n)}$ be the order statistics for the distance values (i.e., sorted into increasing order) we can use a standard probability result for the distribution of order statistics for absolutely continuous random variables, to get the "$p_2$" probabilities for the observed ordered-distances, which are:$^\dagger$

$$\begin{aligned} p_{2}(d_{(i)}) &\equiv \mathbb{P}(D_{(i)} \geqslant d_{(i)}) \\[6pt] &= 1 - \sum_{j=1}^{i-1} {n \choose j} (1 - p_{1}(d_i))^{j} \cdot p_{1}(d_i)^{n-j}. \\[6pt] \end{aligned}$$

This appears to be the formula you are inquiring about. It is derived from the probability rule for order statistics of absolutely continuous random variables.

---

$^\dagger$ Note also that the corresponding marginal densities of the ordered-squared-distances are:

$$\begin{aligned} f_{D_{(i)}}(d_{(i)}) &= \frac{n!}{(i-1)! (n-i)!} \cdot F_D(d_{(i)})^{i-1} \cdot (1-F_D(d_{(i)}))^{n-i} \cdot f_d(d_{(i)}) \\[6pt] &= \frac{n!}{(i-1)! (n-i)!} \cdot (1 - p_{1}(d_i))^{i-1} \cdot p_{1}(d_i)^{n-i} \cdot \text{ChiSq}(d_{(i)}|m). \\[6pt] \end{aligned}$$