Sfakianakis-Verginis quantile estimator

There are dozens of different ways to estimate quantiles. One of these ways is to use the Sfakianakis-Verginis quantile estimator. To be more specific, it’s a family of three estimators. If we want to estimate the $p^\textrm{th}$ quantile of sample $X$, we can use one of the following equations:

where $B_i = B(i; n, p)$ is probability mass function of the binomial distribution $B(n, p)$, $X_{(i)}$ are order statistics of sample $X$.

In this post, I derive these equations following the paper “A new family of nonparametric quantile estimators” (2008) by Michael E. Sfakianakis and Dimitris G. Verginis. Also, I add some additional explanations, reconstruct missing steps, simplify the final equations, and provide reference implementations in C# and R.

Preparation

Consider sample $X = \{ X_1, X_2, \ldots, X_n \}$ ($n \geq 3$). Let $X_{(1)}, X_{(2)}, \ldots, X_{(n)}$ be the order statistics of this sample ($X_{(k)}$ is the $k^\textrm{th}$ smallest element). Now let’s build the following intervals:

where $L$ and $U$ are lower and upper bounds for $X$ values ($L$ and $U$ could be equal to $-\infty$ and $\infty$).

We want to estimate $p^\textrm{th}$ quantile $Q_p$. Obviously, $Q_p$ should belong to one of the $S_i$ intervals (because they cover all possible $Q_p$ values). With the help of $Q_p$, we can introduce variables $\delta_i$:

Since $Q_p$ is the $p^\textrm{th}$ quantile, the probability of $(X_i \leq Q_p)$ is $p$ and the probability of $(X_i > Q_p)$ is $1-p$. Thus, $\delta_i$ belongs to the bernoulli distribution: $\delta_i \sim \textrm{Bernoulli}(p)$.

Next, consider the sum of $\delta_i$ values:

Since $\delta_i$ are independent variables from the Bernoulli distribution, their sum belongs to the binomial distribution: $N \sim \textrm{Binomial(n, p)}$. Now we can get the probability of $Q_p \in S_i$:

Let $Q'_{p,i}$ be a point estimator of $Q_p$ conditioned on the event $Q_p \in S_i$. With the help of $Q'_{p,i}$, we could introduce the Sfakianakis-Verginis quantile estimator:

There are three variations of the Sfakianakis-Verginis quantile estimator depending on the definition of $Q'_{p,i}$.

SV1

The first approach assumes that the $p^\textrm{th}$ quantile is in the middle of the $i^\textrm{th}$ segment:

If we don’t know values of $L$ and $U$, we need a way to define $Q'_{p,i}$ for $i=0$ and $i=n$. The author suggests using these assumptions:

Thus, we have:

Now we can start deriving $\operatorname{SV1}_p$:

It gives us the following expression:

By regrouping the equation members, we get the final result from the paper:

SV2

The second approach assumes that the $p^\textrm{th}$ quantile is in the right endpoint of the $i^\textrm{th}$ segment:

To get the value of $Q'_{p,n}$, we use an assumption from the $\operatorname{SV1}_p$ section:

which can be transformed to

Using $\operatorname{SV2}_p = \sum_{i=0}^n B(i;n,p) Q'_{p,i}$, we get:

SV3

The third approach assumes that the $p^\textrm{th}$ quantile is in the left endpoint of the $i^\textrm{th}$ segment:

To get the value of $Q'_{p,0}$, we use an assumption from the $\operatorname{SV1}_p$ section:

which can be transformed to

Using $\operatorname{SV3}_p = \sum_{i=0}^n B(i;n,p) Q'_{p,i}$, we get:

Simplification

In my opinion, the above equations (that match the original paper) are too bulky. Let’s simplify them! We can notice that we work with the same Binomial distribution $B(n, p)$. What the point of writing the arguments $n$ and $p$ each time? This just makes it difficult to read the equation! (If you read the original paper, you may notice that the authors are also confused with this notation; sometimes they write $B(i;p,n)$ or $B(i;p)$ instead of $B(i;n,p)$.) Let’s denote $B(i;p,n)$ as $B_i$:

Now we could rewrite the first equation in a shorter form:

Currently, the sum part $\sum_{i=2}^{n-1} \frac{B_i+B_{i-1}}{2} X_{(i)}$ enumerates $i$ from $2$ to $n-1$. The part could “absorb” the $\frac{B_i+B_{i-1}}{2} X_{(i)}$ expression for $i=1$ and for $i=n$ from other members of the equation:

In the second equation, we could adjust indexes and enumerate $i$ in range $1..n$ instead of $0..(n-1)$ (to make it consistent with the first and the third equations):

In the third equation, we should just replace $B(i;p,n)$ by $B_i$:

Reference implementation

If you use R, here is the function that you can use in your scripts:

sv1quantile <- function(x, probs) {
  n <- length(x)
  if (n <= 2)
    return(quantile(x, probs))
  x <- sort(x)
  sapply(probs, function(p) {
    B <- function(x) dbinom(x, n, p)
    B(0) * (x[1] + x[2] - x[3]) / 2 +
      sum(sapply(1:n, function(i) (B(i) + B(i - 1)) * x[i] / 2)) +
      B(n) * (-x[n-2] + x[n-1] + x[n]) / 2
  })
}
sv2quantile <- function(x, probs) {
  n <- length(x)
  if (n <= 2)
    return(quantile(x, probs))
  x <- sort(x)
  sapply(probs, function(p) {
    B <- function(x) dbinom(x, n, p)
    sum(sapply(1:n, function(i) B(i - 1) * x[i])) + B(n) * (2 * x[n] - x[n - 1])
  })
}
sv3quantile <- function(x, probs) {
  n <- length(x)
  if (n <= 2)
    return(quantile(x, probs))
  x <- sort(x)
  sapply(probs, function(p) {
    B <- function(x) dbinom(x, n, p)
    sum(sapply(1:n, function(i) B(i) * x[i])) +
      B(0) * (2 * x[1] - x[2])
  })
}

If you use C#, you can take an implementation from the latest nightly version (0.3.0-nightly.90+) of Perfolizer (you need SfakianakisVerginis1QuantileEstimator, SfakianakisVerginis2QuantileEstimator, SfakianakisVerginis3QuantileEstimator).

Conclusion

The Sfakianakis-Verginis quantile estimator is an interesting way to estimate quantiles. In this post, we briefly described the idea and looked at the estimator equations. In future posts, we compare the efficiency of the suggested approach with other quantile estimators.

References

• [Sfakianakis2008]
  Sfakianakis, Michael E., and Dimitris G. Verginis. “A new family of nonparametric quantile estimators.” Communications in Statistics—Simulation and Computation® 37, no. 2 (2008): 337-345.
  https://doi.org/10.1080/03610910701790491