3
$\begingroup$

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?

$\endgroup$

1 Answer 1

2
$\begingroup$

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")

Plot of CDF

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

$\endgroup$
2
  • $\begingroup$ 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)$. $\endgroup$
    – Silverfish
    Commented Feb 12, 2015 at 12:02
  • $\begingroup$ 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)) $\endgroup$
    – Silverfish
    Commented Feb 12, 2015 at 12:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.