How to produce a vector based on two probability distributions?

Question

I am working on distribution on plants and I need to enter the chances that a plant starts reproducing at a certain age in a vector. Only this has to be done for more generations.

Now, Imagine I have a plant that starts reproducing at least at 80 years and at last at 81 years. I assume that the chance for a plant to reach its reproducible age is equally distributed, so the chance a plant starts reproducing at 80 is 0.5 and the chance a plant starts reproducing at 81 is also 0.5. However, I want to know these chances for 2 generations of a plant. So I make a tree diagram, where the chance that 2 generations last 160 years (2*80) is 0.25 (0.5*0.5). The chance that 2 generations last 161 years (80+81 or 81+80) is 0.5 (2*(0.5*0.5)). The chance that 2 generations last 162 years is also 0.25 (0.5*0.5). The vector that I need here is a cumulative vector of chances: c(0.25, 0.75, 1)

However I have a plant that starts reproducing at least at 80 years and at last at 149 years. I assume again that the chances for a plant to reach its reproducible age is equally distributed, so that is 1/70 a plant starts reproducing at 80, 1/70 at 81 etc. I want to know these chances for 20 generations. (You can of course make a tree diagram of this as well, but that will be one enormous tree diagram). This will become a vector of 1400 values, where the first value represents the chance that all generations start reproducing young and that is in (80*20)= 1600 years and the last value represents the chance that all genereations start reproducing at ist oldest age (150*20)= 3000 years. In this vector each value stands for the chance in one year (3000-1600=1400).

So here's my question: Is there a way in R to calculate such a cumulative vector of chances?

Answer 1

The successive generations correspond to adding random variables representing the time of reproduction. Adding variables is a computational convolution. It is most efficiently done with a Fast Fourier Transform. R carries this out with convolve. The resulting probabilities can be added with cumsum.

cdf <- function(x, n=2) {
  if (n < 2) return(cumsum(x))
  y0 <- rev(x); y <- x
  for (i in 1:(n-1)) y <- convolve(y, y0, type="open")
  cumsum(y)
}

The first argument is a vector of probabilities $x = (x_0, x_1, \ldots)$. The coefficients $x_i$ are the chances of reproducing at year $i$. The second argument is the number of generations, $n$, defaulting to $2$. Let's apply this solution to the example in the question:

plant <- c(rep(0,80), rep(1/70, 70))
plot(cdf(plant, 20), type="l", xlab="Year")

The calculation takes about $0.04$ seconds. It is still inefficient: direct use of fft ought to speed it up by an order of magnitude. But that requires additional manipulation and testing, so for reliability, convolve is a good choice for small problems like this one.

---

Incidentally, after 20 generations the distribution will become extremely close to normal. The original one is uniform over integers from $80$ through $149$, whence it has a mean of $114.5$, and will have a variance extremely close to that of a truly uniform distribution of width $70$, equal to $70^2/12$. Indeed, if you overplot the sum of 20 copies of this normal distribution you will see essentially no difference:

curve(pnorm(x, 20*(114.5), 70*sqrt(20/12)), add=TRUE)

Whence another way to obtain the desired probabilities is

z <- pnorm(seq(from=1/2, along.with=y), 20*(114.5), 70*sqrt(20/12))

These values differ from the output of cdf in this example by no more than $0.0014$.