# Inter-arrival time of subsampled Poisson point process

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!

• Since $n$ is random, if $m$ is constant there will be a non-zero probability that $n<m$, in which case you can not keep $m$ points.
– Yves
Sep 8, 2015 at 15:15
• $n$ is fixed, not random. I should emphasize that I mean a Poisson point process. I will have $n$ points in the end, it's the time it'll take me that is random.
– lum
Sep 8, 2015 at 15:25
• This could be told for clarity: so $t$ is the waiting time for the $n$-th arrival of the PP. Both $n$ and $m <n$ are fixed.I guess that the choice of the $m$ points among the $n$ is made independently of the process.
– Yves
Sep 8, 2015 at 15:46
• A 'new' interarrival (IE) is the sum of $K$ 'old' IE where $K$ is a r.v. with integer values from $1$ to $n-m+1$, which is independent from the old IE. So the new IE is a mixture of gammas with shape $k$ and rate $\lambda$. There remains to compute the weights $\text{Pr}\{K = k\}$ which may relate to the runs distribution. I think that $E(K) = (n-1) / (m-1)$ because $m-1$ new IE are made from the $n-1$ old ones.
– Yves
Sep 9, 2015 at 6:40
• There are errors in my comment, see my detailed answer.
– Yves
Sep 9, 2015 at 15:44

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)

• Yves, thank you for your very thorough answer! I will have a close look as soon possible.
– lum
Sep 10, 2015 at 7:32
• Yves, a much belated feedback: your answer checks out! Once again, thank you very much for the help. Minor typo: on line 4, it should be $I_{m+1} := n+1$ (and not $I_m$.)
– lum
Oct 21, 2015 at 11:41