# How to produce a vector based on two probability distributions?

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?

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$.

• I'll leave a quick note about the convolution code for later readers, as this Q&A is the top google hit for "convolve probability in R". For the common example of the distribution of the sum of two independent dice, one can find the non-zero probabilities by cdf(rep(1/6, 6), 2) but NB the first element is $F_{\sum{X_i}}(2)$ (as is clear from the answer; you might want to multiply the ouput by 36 to see this more clearly). Similarly the first element of cdf(rep(1/6, 6), 8) is $F_{\sum{X_i}}(7)$. Commented Feb 12, 2015 at 12:02
• Setting the first element of x to $P(X_i=0)$, as the answer suggests, includes the zero probabilities as well. In this case the $i^{th}$ element of the output of cdf equals $F(i-1)$, so the first element is $F(0)$ and so on. The output from cdf(c(0,rep(1/6, 6)), 2) starts [1] 9.474259e-18 5.718124e-18 2.777778e-02 8.333333e-02 1.666667e-01 yet we know from theory that $F(0)=F(1)=0$ (the minimum score from rolling two dice is two). This slight discrepancy comes from the FFT calculations; the R help page on convolve suggests using zapsmall, e.g. zapsmall(cdf(c(0,rep(1/6, 6)), 2)) Commented Feb 12, 2015 at 12:11