6
$\begingroup$

[Edit: for posterity, I'm adding two small comments to the code explaining how to fix it, in light of Iosef Pinelis' answer below. Look for "Should be:" to find the corrections.]

Suppose we draw a sample of $n$ i.i.d. observations from some continuous distribution on $\mathbb{R}$. Sort these observations from smallest to largest and call them $x_1,\ldots,x_n$. We consider the quantile corresponding to each observation, $q_1,\ldots,q_n$. We are interested in the probability that the $q_i$ are all simultaneously close to $i/n$.

The DKW inequality shows that $$ \Pr\left[\max_{i=1,\ldots,n}\left|q_i-\frac{i}{n}\right|>\varepsilon\right] \leq 2 e^{-2n\varepsilon^2} $$ Note that this result holds for all $n$ and there are no hidden constants.

It is extremely straightforward to run a Monte Carlo simulation of this process. The uniform distribution on $[0,1]$ is particularly convenient to use since each obersvation in the sample equals its own quantile (i.e., $x_i=q_i$). I've shared some Python code; 1 million Monte Carlo trials take about 10 seconds.

My simulation appears to violate the DKW inequality. Question: What's going on?

# Double-check the DKW Inequality
import numpy as np

################
# Set parameters
# 
# "prob_target" is the probability of *satisfying* the bound
# (i.e., that the gap is less than epsilon)
prob_target = 0.99
# We consider the CDF of "n" samples
n=20
# We are going to use "num_trials" Monte Carlo runs, where each run
# consists of "n" draws from a Uniform distribution
num_trials=1000000
print(f"Parameters: n={n}, # trials={num_trials}, target probability={prob_target}")

#######################
# Compute DKW threshold
epsilon_DKW = np.sqrt(np.log(2.0 / prob_target) / (2.0 * float(n)))
#   Should be:
# epsilon_DKW = np.sqrt(np.log(2.0 / (1-prob_target)) / (2.0 * float(n)))
# Double-check that we actually computed the intended probability;
# hopefully prob_target == prob_DKW
prob_DKW = 2 * np.exp(-2 * n * (epsilon_DKW**2))
#   Should be:
# prob_DKW = 1 - 2 * np.exp(-2 * n * (epsilon_DKW**2))
assert np.isclose(prob_target, prob_DKW)
print(f"Computed parameters: epsilon_DKW ={epsilon_DKW:.6f}, prob_DKW={prob_DKW:.6f}")

########################
# Monte Carlo simulation
disparity_list = np.zeros(num_trials)
ecdf = np.arange(1, n+1)/n  # = [1/n, 2/n, ..., n/n]
for trial in range(num_trials):
    data = np.random.uniform(size=n)
    data = np.sort(data)
    quantiles = data  # ...because uniform distribution on [0,1]
    worst_disparity = np.max(np.abs(quantiles-ecdf))

    disparity_list[trial] = worst_disparity
# Compute fraction of trials with disparity below the DKW bound
prob_true = np.sum(disparity_list < epsilon_DKW) / num_trials
# Compute the *actual* bound such that prob_DKW fraction of
# trials have disparity below the bound
epsilon_true = np.quantile(disparity_list, prob_DKW)

###############
# Print results
print(f"Measured threshold:  epsilon_best={epsilon_true:.6f}")
if prob_true < prob_DKW:
    print("\nWe have a problem.")
    print(f"DKW promises success probability at least {prob_DKW:.6f}, ")
    print(f"but we only observe probability           {prob_true:.6f}")
    print(f"The required epsilon is {epsilon_true / epsilon_DKW:.6f}x larger than DKW!")
else:
    print(f"DKW Inequality works!  Success probability={prob_true:.6f}>={prob_DKW:.6f}")

Here's the result of the run:

Parameters: n=20, # trials=1000000, target probability=0.99
Computed parameters: epsilon_DKW =0.132589, prob_DKW=0.990000
Measured threshold:  epsilon_best=0.335025

We have a problem.
DKW promises success probability at least 0.990000, 
but we only observe probability           0.333689
The required epsilon is 2.526787x larger than DKW!
$\endgroup$
7
  • $\begingroup$ In your definition of prob_true did you have a reversed inequality sign? $\endgroup$
    – John Jiang
    Commented Dec 31, 2022 at 6:21
  • $\begingroup$ Ah, @JohnJiang, I believe that it's correct as it stands, but that touches on the mistake I made. (See the comments below Iosef Pinelis' solution below.) But the ambiguity there is what killed me; I'll add a comment to make it clearer. $\endgroup$ Commented Dec 31, 2022 at 14:32
  • $\begingroup$ The standard way in which the word "sample" is used in statistics would refer to "a sample of $n$ i.i.d. observations" rather than to "$n$ i.i.d. samples". For example, the thing that in standard usage is called a "two-sample t-test" obviously does not rely on just two independent observations, and what would the term "sample size" mean? $\endgroup$ Commented Dec 31, 2022 at 16:07
  • $\begingroup$ Thank you, @MichaelHardy, I've fixed that. $\endgroup$ Commented Dec 31, 2022 at 16:09
  • $\begingroup$ $\frac in$ looks biased upwards. $q_i \sim \text{Beta}(i,n-i+1)$ so the expected value of $q_i$ is $\frac{i}{n+1}$, the most-likely/mode/highest-density value is $\frac{i-1}{n-1}$ for $n>1$, and the median is about $\frac{i-\frac13}{n+\frac13}$ for $1 < i <n$ $\endgroup$
    – Henry
    Commented Jan 3, 2023 at 0:05

1 Answer 1

5
$\begingroup$

The DKW inequality can be stated as follows: \begin{equation} p_n(z):=P(\|S_n\|_\infty>z)\le 2e^{-2z^2} \tag{1}\label{1} \end{equation} for $z>0$, where \begin{equation} S_n(t):=\frac1{\sqrt n}\,\sum_{i=1}^n Y_i(t),\quad Y_i(t):=1(X_i\le t)-t, \end{equation} and the $X_i$'s are iid random variables (r.v.'s) each uniformly distributed on $[0,1]$.

You seem to say that your simulations appear to be suggesting that the constant factor $2$ in the exponent of $e^{-2z^2}$ in the bound in \eqref{1} must be replaced a constant factor that is about $2.5$ times worse than the constant factor $2$.


This suggestion is extremely unlikely to be true, given the very long and rich history of the DKW inequality. Apparently, this history begins with Kolmogorov and Smirnov, who showed that
\begin{equation} p_n(z)\to p(z):=2\sum_{k=1}^\infty(-1)^{k-1}e^{-2k^2z^2} \end{equation} for each $z>0$ as $n\to\infty$. Note that $p(z)\le 2e^{-2z^2}$ for all $z>0$ and $p(z)\sim 2e^{-2z^2}$ as $z\to\infty$. This strongly suggests that $2$ is the right (and the best possible) constant factor $2$ in the exponent of $e^{-2z^2}$ in \eqref{1}.

One may also note here the iid random functions $Y_i$ are zero-mean, each with the covariance function $[0,1]^2\ni(s,t)\mapsto EY_i(s)Y_i(t)=\min(s,t)-st$. So, the mean and covariance functions of each random function $Y_i$ are the same as those of the Brownian bridge. By an appropriate functional central limit theorem, it is at least plausible that the random function $S_n$ converges in an appropriate sense to a Gaussian process. Since the distribution of a Gaussian process is determined by its mean and covariance functions, it follows that the random function $S_n$ converges in an appropriate sense to a Brownian bridge $B^0$. So, it is likely that \begin{equation} p_n(z)\to P(\|B^0\|_\infty>z)=p(z); \tag{2}\label{2} \end{equation} concerning the latter equality, cf. e.g. this or this. Of course, \eqref{2} is, not only plausible, but true and well known. In particular, in view of (say) [Talangrand's Theorem 1.1] with $\mathcal F=\{1_{[0,t]}\colon t\in[0,1]\}$, $\|S_n-B^0\|_\infty$ converges to $0$ in probability as $n\to\infty$.

So, again, it is extremely unlikely that the DKW inequality is false. It is much more likely that the problem is with the software you were using and/or with your code, which I cannot read.


Below is an image of a Mathematica notebook confirming your inequality:

enter image description here

$\endgroup$
7
  • $\begingroup$ You could also write the function as (RandomReal[{0,1},n] // Sort - Range[n] / n) // Abs // Max $\endgroup$
    – user44143
    Commented Dec 30, 2022 at 2:25
  • $\begingroup$ @MattF. : Thank you for your comment. I have used Mathematica for many years, but have never become good at using it. $\endgroup$ Commented Dec 30, 2022 at 4:48
  • $\begingroup$ I have seen you get plenty of good results quickly with Mathematica! For me, the key to using it well is knowing that you can write most loops more elegantly with built-in functions. Here Range is helpful, but I also frequently find second-order functions like Fold/FoldList, Nest/NestList, Map/MapAt, Outer, Thread to be useful, even if figuring out which one I need requires a pause to review some examples first. $\endgroup$
    – user44143
    Commented Dec 30, 2022 at 14:36
  • $\begingroup$ @MattF. : I have used these commands, maybe a couple times each (except for Map, which I have used quite a few times). Yet, I keep forgetting how to use those commands or even their names. $\endgroup$ Commented Dec 30, 2022 at 14:48
  • $\begingroup$ Thank you so much for this thoughtful answer. In case it wasn't clear from my question, obviously I assume that the DKW inequality is correct and something about my own approach was in error. I will now try to deduce enough Mathematica to read your code, compare it with my own, and look for discrepancies. $\endgroup$ Commented Dec 30, 2022 at 14:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.