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

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

• It's not clear if you want the number of individuals moving to be exactly 5 or not. If you do want exactly 5, then couldn't you just choose 5 randomly and then randomize their move times however you want? – soakley Jan 11 at 17:47
• @soakley no, it should be probability-based, so at the end of the time period, approximately 95 individuals remain – jayb Jan 12 at 11:45
• I would still divide the processes. Simulate 100 realizations from a Bernoulli random variable with success probability 0.05 - that gives you who moves away and will average 5 people. Then, if you need a time when they move away (I don't know if you need that or not), you can simulate that for those that have been selected as movers. – soakley Jan 12 at 13:57

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.

• Thanks for the detailed answer. Is it not possible to estimate a daily probability of moving for each individual then? Maybe I've misunderstood, but it seems like your method estimates the time when d individuals move, not the probability that any individual will move on a given day – jayb Jan 12 at 16:39
• This analysis produces the answers to all such questions of probabilities. Indeed, at the outset it states the probability per day in terms of the rate, which is a function of how many people are left. There is no estimation involved: given your assumptions, the calculations are exact. Now, some questions do require a calculation based on this formulation. For instance, if you want to compute an unconditional probability per day, then you will need to compute that integral. But as far as I can tell, that's not what you asked: your question asks about simulation. – whuber Jan 12 at 18:11
• Thanks for clarifying. It is an unconditional probability per day that I'm after. Could you expand your response to include this? – jayb Jan 12 at 19:52