R - generate deviate using probability based on known occurrence of event

Question

In R, I want to simulate an event using a probability $p$ based on a known occurrence $k$ of the event.

If I have a population of $n = 100$ individuals in a given area, and I know over a period of time $t$ (e.g. 365-days) that 5 of these individuals will move away from that area, what is the best way to simulate this?

Each day, for each individual, a binomial deviate is drawn indicating whether or not the individual will move on that day. An individual can only move away from the area once, and once moved cannot move back - therefore, the population will decrease over time. At the end of the time period (e.g. 365-days), there should be approximately 95 individuals remaining.

At first, I thought the daily probability* of an event such as this should be derived as $ p = 1 - (1 - \frac{k}{n})^{1/t} $, where $k = 5, n = 100, t = 365$, and simulated in R using a binomial draw rbinom() with $p$ probability. But I'm not sure that this is best way and would appreciate some guidance.

Edit *unconditional probability - as clarified by @whuber

Answer 1

Assume the actions of these individuals are independent. (If you don't, you need to specify--mathematically--how they are interdependent.)

Assume that the chance any individual leaves during any time period does not change over time. (This can be modified, but then you will have to specify how these chances vary over time.)

These assumptions imply that the time $T_i$ when individual $i$ leaves is a random variable with an exponential distribution. You can make progress towards a solution by supposing every individual has the same propensity to leave: that is, the rate parameters of these exponential distributions are all the same. (If you don't want to assume that, you must make some assumption about the distribution of rate parameters in the population.) Let this common rate be $\kappa.$

These assumptions imply that when there are $n$ people, the time until the next person leaves (the smallest of the $T_i$) has an exponential distribution with rate $n\kappa.$ Its expectation is $1/(n\kappa).$ After that $n-1$ people remain, so the time until the second person leaves has an exponential distribution with rate $(n-1)\kappa.$ And so on. Consequently, the distribution of the time taken until the fifth person leaves is the sum of (independent) exponential distributions of rates $n\kappa,$ $(n-1)\kappa,$ and on on down to $(n-4)\kappa.$

The expectation of the sum is the sum of the expectations of the terms, equal to

$$\mu = \frac{1}{n\kappa} + \frac{1}{(n-1)\kappa} + \cdots + \frac{1}{(n-4)\kappa}.$$

Solving for $\kappa$ gives

$$\kappa = \frac{H(n) - H(n-5)}{\mu}$$

where $H(n)$ is the Harmonic Number.

In your problem you stipulate $\mu=5$ people per year and $n=100.$ The corresponding value of $\kappa$ is $0.01020621.$

In general (depending on $n$ and $\mu$) two simulation approaches look feasible. One is to simulate the sequence of exponential variates and add them. The other is to compute the distribution once and for all. According to the solution at https://stats.stackexchange.com/a/72486/919, the distribution (of the time until the fifth person leaves) is a mixture of exponential distributions. This isn't any simpler or more efficient to compute, so let's go with the first method.

In R, for instance, the time until d people out of n will leave can be simulated in two lines: one to compute $\kappa$ and the other to generate the exponential variates and sum them:

kappa <- sum(1/seq(n, n-d+1, by=-1))
x <- colSums(matrix(rexp(1e6 * d, rate=seq(n, n-d+1, by=-1)*kappa), d))

(The cumulative column sums in this matrix would indicate the specific times at which the first, second, ..., through the $d^\text{th}$ person leaves.)

This example generates a dataset of 1e6 values. Here is their histogram:

hist(x, freq=FALSE, breaks=200, xlim=c(0, 3.5), border="Gray",
     main=paste("Histogram of Time for", d, "People to Leave"))

The red curve was computed using the partial fractions formula previously quoted. Clearly it agrees with the simulation, confirming the correctness of both approaches.