Inter-arrival time of subsampled Poisson point process

Question

Suppose that I draw $n$ points from a Poisson point process of rate $\lambda$, i.e. with inter-arrival times distributed i.i.d $\sim \text{Exp}(\lambda)$.

Now suppose that I choose $m < n$ of these points uniformly at random, and throw away the rest.

My question is as follows: What can be said about the inter-arrival times in this modified process? I.e. about the inter-arrival time between the $m$ remaining points?

A couple of related facts that I know:

1. if each point in the original point process is kept independently with probability $p$, then the "thinned" process has inter-arrival times i.i.d $\sim \text{Exp}(p \lambda)$
2. If the $m$ points are chosen i.i.d with a uniform distribution over some interval $[0,t]$, and then ordered, the inter-arrival times would be Beta distributed.

Somehow, neither of these observations seems to fit my case. Any pointer is appreciated!

Answer 1

Let $n$ and $m<n$ be fixed. Consider first the random sampling without replacement of $m$ values in $\{1,\,2,\, \dots,\,n\}$, leading to $m$ order statistics $I_1 < I_2 < \dots < I_m$ for the retained numbers and $m-1$ spacings $K_j:= I_{j+1} -I_{j}$ for $j=1$ to $m-1$. With the conventions $I_0:=0$ and $I_{m+1}:=n+1$, we define as well $K_0=I_1$ and $K_m=n+1-I_m$. It is not too difficult to see that the $m + 1$ random variables $K_j$ have the same distribution, namely the distribution of $I_1$. Moreover $$ \text{Pr}\{ I_1 = k \} = \frac{{n-k}\choose{m-1}}{{n}\choose{m} } := p_k \qquad (1 \leq k \leq n-m +1), $$ since the number of 'favourable outcomes' is the number of choices of $m-1$ items among $n-k$. This allows the computation of the $p_k$ by the recursion $$ p_{k+1} = \frac{ n- k -m +1}{n-k} \, p_{k} $$ for $k = 1$ to $n-m$, starting from a positive $p_1$ and then normalising to have total mass unity. Note that $\mathbb{E}[K] = (n+1)/(m+1)$.

Now let us come back to the Poisson process and the "ensemble thinning" of the question. If the $n$ arrivals are at times $X_1 < X_2 < \dots < X_n$ the $m$ selected arrivals are $Y_j:=X_{I_j}$ where the $I_j$ are as before. The interarrival $W_j:=Y_{j+1}-Y_{j}$ is the sum of the random number $K_j$ of independent exponential interarrivals $X_{i+1}-X_{i}$. Conditional on $K_j=k$ the r.v. $W_j$ is Gamma with shape $k$ and rate $\lambda$. The unconditional distribution of a $W$ is a mixture $$ W \sim \sum_{k=1}^{n-m+1} p_k \, \text{Gam}(k,\,\lambda). $$ We implicitely supposed that the choice of the $m$ retained arrivals is done independently of the $X_i$. Note that the $W_j$ are not independent because the $K_j$ are not independent.

n <- 30; m <- 12 

## Monte-Carlo simulation with 'N' replications ('lambda' is 1) 
N <- 30000 
set.seed(1234)

## compute the sample distributions for the new interarrivals 'W'
## and the intersampled 'K'
W <- matrix(NA, nrow = N, ncol = m - 1)
K <- matrix(NA, nrow = N, ncol = m + 1)
colnames(K) <- paste("K", 0:m, sep = "_")
colnames(W) <- paste("W", 1:(m - 1), sep = "_")

for (samp in 1:N) {
    X <- cumsum(rexp(n))                    ## old arrivals
    Isamp <- sort(sample(1:n, size = m))    ## indices for new arrivals
    K[samp, ] <- diff(c(0, Isamp, n + 1))   ## compute 'K' for checks
    W[samp, ] <- diff(X[Isamp])             ## new interarrivals
}

## compute the exact distribution of the r.vs 'K' 
p <- rep(NA, n - m + 1)
p[1] <- 1
for (k in 1:(n - m)) {
    p[k + 1] <- p[k] * (n - k - m + 1) / (n - k)
}
p <- p / sum(p)

## compute the exact expectation of 'K' for check
EK <- sum(p * seq(from = 1, to = n - m + 1, by = 1))

## compare distribution of W_j, 'j' is between 0 and 'm'
j  <- 8
tab <- table(factor(K[ , j + 1], levels = 1:(n - m + 1)))
compar <- cbind(Sim = tab, Exact = round(p * N)) 
cols <-  c("SteelBlue", "orangered")
barplot(t(compar), beside = TRUE,
        legend.text = c(sprintf("Simul. K_%d", j), "Exact"),
        main = sprintf("Mixture weights p_k for n = %d, m = %d", n, m),
        col = cols)

## compare distributions functions for interarrivals
F_Sim <- ecdf(W[ , 1])

qx <- quantile(F_Sim)
x <- seq(from = qx["0%"], to = qx["100%"], length.out = 300)
F_Exact <- rep(0, length(x))
for (k in 1:(n - m + 1)) {
    F_Exact<- F_Exact+ p[k] * pgamma(x, shape = k)
}


plot(F_Sim, col = cols[1],
     main = "distribution function of the 'new' interarrival W_1",
     xlab = "w", ylab = "Fn(w)", lwd = 2)
lines(x, F_Exact, col = cols[2], lty = "dashed", lwd = 2)
legend("center", legend = c("Simul. (ECDF)", "Exact (mixt. of gammas)"),
       lty = c("solid", "dashed"), lwd = 2.5, col = cols)