17
$\begingroup$

This problem is related to my lab's research in robotic coverage:

Randomly draw $n$ numbers from the set $\{1,2,\ldots,m\}$ without replacement, and sort the numbers in ascending order. $1\le n\le m$.

From this sorted list of numbers $\{a_{(1)},a_{(2)},…,a_{(n)}\}$, generate the difference between consecutive numbers and the boundaries: $g = \{a_{(1)},a_{(2)}−a_{(1)},\ldots,a_{(n)}−a_{(n-1)},m+1-a_{(n)}\}$. This gives $n+1$ gaps.

What is the distribution of the maximum gap?

$P(\max(g) = k) = P(k;m,n) = ? $

This can be framed using order statistics: $P(g_{(n+1)} = k) = P(k;m,n) = ? $

See link for the distribution of gaps, but this question asks the distribution of the maximum gap.

I'd be satisfied with the average value, $\mathbb{E}[g_{(n+1)}]$.

If $n=m$ all the gaps are size 1. If $n+1 = m$ there is one gap of size $2$, and $n+1$ possible locations. The maximum gap size is $m-n+1$, and this gap can be placed before or after any of the $n$ numbers, for a total of $n+1$ possible positions. The smallest maximum gap size is $\lceil\frac{m-n}{n+1}\rceil$. Define the probability of any given combination $T= {m \choose n}^{-1}$.

I've partially solved the probability mass function as $ P(g_{(n+1)} = k) = P(k;m,n) = \begin{cases} 0 & k < \lceil\frac{m-n}{n+1}\rceil\\ 1 & k = \frac{m-n}{n+1} \\ 1 & k = 1 \text{ (occurs when $m=n$)} \\ T(n+1)& k = 2 \text{ (occurs when $m=n+1$)} \\ T(n+1)& k = \frac{m-(n-1)}{n} \\ ? & \frac{m-(n-1)}{n} \le k \le m-n+1 \\ T(n+1)& k = m-n+1\\ 0 & k > m-n+1 \end{cases} \tag{1} $

Current work (1): The equation for the first gap, $a_{(1)}$ is straightforward: $$ P(a_{(1)} = k) = P(k;m,n) = \frac{1}{{m \choose n}} \sum_{k=1}^{m-n+1} {m-k-1 \choose n-1} $$ The expected value has a simple value: $\mathbb{E}[P(a_{(1)})] = \frac{1}{ {m \choose n}} \sum_{k=1}^{m-n+1} {m-k-1 \choose n-1} k = \frac{m-n}{1+n}$. By symmetry, I expect all $n$ gaps to have this distribution. Perhaps the solution could be found by drawing from this distribution $n$ times.

Current work (2): it is easy to run Monte Carlo simulations.

simMaxGap[m_, n_] := Max[Differences[Sort[Join[RandomSample[Range[m], n], {0, m+1}]]]];
m = 1000; n = 1; trials = 100000;
SmoothHistogram[Table[simMaxGap[m, n], {trials}], Filling -> Axis,
Frame -> {True, True, False, False},
FrameLabel -> {"k (Max gap)", "Probability"},
PlotLabel -> StringForm["m=``,n=``,smooth histogram of maximum map for `` trials", m, n, trials]]

enter image description here

$\endgroup$
9
  • 1
    $\begingroup$ With these conditions you must have n<=m. I think you want g={a_(1), a_(2)-a_(1),..., a_(n)-a_(n-1)}. Does randomly select mean selecting each number with probability 1/m on the first draw? Since you do not replace the probability would be 1/(m-1) on the second and so on down to 1 on the mth draw if n=m. If n<m this would stop earlier with the last draw having probability 1/(m-(n-1)) on the nth draw. $\endgroup$ Dec 31, 2016 at 19:34
  • 2
    $\begingroup$ Your original description of $g$ made no sense, because (I believe) you transposed two of the subscripts. Please verify that my edit conforms with your intention: in particular, please confirm that you mean for there to be $n$ gaps, of which $a_{(1)}$ is the first. $\endgroup$
    – whuber
    Dec 31, 2016 at 19:55
  • 1
    $\begingroup$ @gung I think this is research, rather than self-study $\endgroup$
    – Glen_b
    Jan 1, 2017 at 18:20
  • 1
    $\begingroup$ I think your minimum and maximum gap sizes should be $1$ and $m-n+1$. The minimum gap size is when consecutive integers are chosen, and the maximum gap size occurs when you select $m$ and $n-1$ first integers $1,\dots,n-1$ (or $1$ and $m-n+2,\dots,m$) $\endgroup$ Jan 1, 2017 at 23:12
  • 1
    $\begingroup$ Thank you Michael Chernick and probabilityislogic, your corrections have been made. Thank you @whuber for making the correction! $\endgroup$ Jan 2, 2017 at 13:32

2 Answers 2

10
$\begingroup$

Let $f(g;n,m)$ be the chance that the minimum, $a_{(1)}$, equals $g$; that is, the sample consists of $g$ and an $n-1$-subset of $\{g+1,g+2,\ldots,m\}$. There are $\binom{m-g}{n-1}$ such subsets out of the $\binom{m}{n}$ equally likely subsets, whence

$$\Pr(a_{(1)}=g = f(g;n,m) = \frac{\binom{m-g}{n-1}}{\binom{m}{n}}.$$

Adding $f(k;n,m)$ for all possible values of $k$ greater than $g$ yields the survival function

$$\Pr(a_{(1)} \gt g) = Q(g;n,m)= \frac{(m-g)\binom{m-g-1}{n-1}}{n \binom{m}{n}}.$$

Let $G_{n,m}$ be the random variable given by the largest gap:

$$G_{n,m} = \max\left(a_{(1)}, a_{(2)}-a_{(1)}, \ldots, a_{(n)}-a_{(n-1)}\right).$$

(This responds to the question as originally framed, before it was modified to include a gap between $a_{(n)}$ and $m$.) We will compute its survival function $$P(g;n,m)=\Pr(G_{n,m}\gt g),$$ from which the entire distribution of $G_{n,m}$ is readily derived. The method is a dynamic program beginning with $n=1$, for which it is obvious that

$$P(g;1,m) = \Pr(G_{1,m} \gt 1) = \frac{m-g}{m},\ g=0, 1, \ldots, m.\tag{1}$$

For larger $n\gt 1$, note that the event $G_{n,m}\gt g$ is the disjoint union of the event

$$a_{1} \gt g,$$

for which the very first gap exceeds $g$, and the $g$ separate events

$$a_{1}=k\text{ and } G_{n-1,m-k} \gt g, \ k=1, 2, \ldots, g$$

for which the first gap equals $k$ and a gap greater than $g$ occurs later in the sample. The Law of Total Probability asserts the probabilities of these events add, whence

$$P(g;n,m) = Q(g;n,m) + \sum_{k=1}^g f(k;n,m) P(g;n-1,m-k).\tag{2}$$

Fixing $g$ and laying out a two-way array indexed by $i=1,2,\ldots,n$ and $j=1,2,\ldots,m$, we may compute $P(g;n,m)$ by using $(1)$ to fill in its first row and $(2)$ to fill in each successive row using $O(gm)$ operations per row. Consequently the table can be completed in $O(gmn)$ operations and all tables for $g=1$ through $g=m-n+1$ can be constructed in $O(m^3n)$ operations.

Figure

These graphs show the survival function $g\to P(g;n,64)$ for $n=1,2,4,8,16,32,64$. As $n$ increases, the graph moves to the left, corresponding to the decreasing chances of large gaps.

Closed formulas for $P(g;n,m)$ can be obtained in many special cases, especially for large $n$, but I have not been able to obtain a closed formula that applies to all $g,n,m$. Good approximations are readily available by replacing this problem with the analogous problem for continuous uniform variables.

Finally, the expectation of $G_{n,m}$ is obtained by summing its survival function starting at $g=0$:

$$\mathbb{E}(G_{n,m}) = \sum_{g=0}^{m-n+1} P(g;n,m).$$

Figure 2: contour plot of expectation

This contour plot of the expectation shows contours at $2, 4, 6, \ldots, 32$, graduating from dark to light.

$\endgroup$
1
  • $\begingroup$ Suggestion: line "Let $G_{n,m}$ be the random variable given by the largest gap:", please add the last gap of $m+1-a_{n}$. Your expectation plot matches my Monte Carlo simulation. $\endgroup$ Jan 2, 2017 at 18:16
1
$\begingroup$

Using the notation in the question, this answer provides a closed-form solution for the number of $(n+1)$-compositions of $m+1$ constrained by $\max(g)\leq k$:

$$\sum_{j=0}^{n+1}(-1)^j\binom{n+1}{j}\binom{m-jk}{n}$$

and $\big\lceil{\frac{m+1}{n+1}}\big\rceil\leq{k}\leq{m-n+1}$

The CDF of $k$ is then

$$\frac{\sum_{j=0}^{n+1}(-1)^j\binom{n+1}{j}\binom{m-jk}{n}}{\binom{m}{n}}$$

The PDF implemented as an R function:

dmaxgap <- function(m, n) {
  w <- ceiling((m + 1)/(n + 1)):(m - n + 1L)
  j <- 0:(n + 1L)
  setNames(diff(c(0L, colSums(rep(c(1, -1), length.out = n + 2L)*choose(n + 1L, j)*choose(pmax(m - outer(j, w), 0L), n))))/choose(m, n), w)
}

Compare the results to $10^6$ simulation replications:

library(Rfast)
m <- 1e3L
n <- 6L

kprobs <- dmaxgap(m, n)

k <- colMaxs(diff(colSort(rbind(0L, replicate(1e6, sample(m, n)), m + 1L))), TRUE)
plot(x <- sort(unique(k)), ecdf(k)(x), type = "l", col = "blue", xlab = "w", ylab = "CDF")
lines(as.integer(names(kprobs)), cumsum(kprobs), col = "orange")

enter image description here

The two lines are nearly indistinguishable.

$\endgroup$

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.