Suppose a female in a given species has, on a given year $x$, PMF $Bernoulli(f(x))$ of finding a mate given she doesn't yet have one. Once she finds a mate, she doesn't need another.

This same female has, on a given year $x$, PMF $Bernoulli(g(x))$ of conceiving and giving birth to 1 offspring, given she already has a mate.

The question is: by age $x$, what is the probability that a female will have at least $k$ offspring?

Is there an analytical solution to this problem when $f,g$ are constants? What is a suitable numerical method for solving this problem when $f,g$ are arbitrary functions?

  • $\begingroup$ I'm missing something. Isn't this just the binomial distribution with chance of reproductive success in a given year equal to f(x)*g(x) for that year? $\endgroup$
    – user1566
    Sep 27 '16 at 14:52
  • $\begingroup$ Not quite. Have a look at Antoni's answer. $\endgroup$
    – Alexandre
    Sep 28 '16 at 19:49

The probability of mating at any given year is $\small \Pr(\text{mate})=m$, and the probability of offspring given a mate has been found is $\small \Pr( \text{single offspring} \vert \text{mate} ) =o$ and it does not change after the mate is found.

The probabilities of getting $k$ offspring after $x$ years depends on the year at which the mate came into the picture.

If we label the beginning of the "experiment" as year zero, $\small \text{Yr}=0$, the probability of having a single offspring at year $0$ is simply $\small \Pr(\text{mate } \cap \text{ offspring})= \Pr(\text{offspring}\vert \text{mate})\Pr(\text{mate})=o\times m.$

At year $1$ (one year later) the probability of a single offspring is simply going to be $\small\Pr(\text{offspring}\vert \text{mate})=o$: the participation of a mate is now guaranteed.

So we just need to treat every year after the appearance of the mate as a binomial, focusing on the year the mate is "realized": the probability that year $1$ is the first year with a mate is calculated based on the probability of absence of mate at year zero, $1-m$, as $\small \Pr_{\text{Yr=1}}(\text{mate})=(1-m)\,m$. At year two, $\small \Pr_{\text{Yr=2}}(\text{mate})=(1-m)^2\,m$; and in general, $\small \Pr_{\text{Yr=yr}}(\text{mate})=(1-m)^{\text{yr}}\,m.$

So we have that the probability of having $k$ offspring by year $x$ and that the mate appears at year $\small \text{Yr=yr}$ is:

\begin{align} \small \Pr(\text{Off= }k\text{ by time } x \cap \text{mate @ Yr=yr})&=\small \Pr(\text{Off= }k\,\vert \,\text{mate @ Yr=yr})\,\Pr(\text{mate @ Yr=yr})\\ &=\small\left({x-\text{yr}+1\choose k}\,o^k\,(1-o)^{x-\text{yr}+1 -k}\right)\, m(1-m)^{\text{yr}}\tag{*} \end{align}

Now, the probability of getting $k$ offspring by year $x$ can happen in any of these scenarios where the mate appears in some year or other, so:

$$\small \Pr(\text{Off}=k\text{ by time } x)=\sum_{\text{Yr=0}}^{x-k+1} m(1-m)^{\text{yr}}{x-\text{yr}+1\choose k}\,o^k\,(1-o)^{x-\text{yr}+1 -k}$$

Notice that the mate cannot appear any later than $x -k + 1$ years if we need $k$ offspring, explaining the upper limit of the sum.

Finally, your question is $\text{at least }k$, inviting another summation:

$$\small\Pr(\text{Off}\geq k\text{ by time } x)= \sum_{\text{Yr=0}}^{x-k+1}\,\, \sum_{\text{K}=k}^{x-\text{yr}+1} \,\, m(1-m)^{\text{yr}}{x-\text{yr}+1\choose k}\,o^k\,(1-o)^{x-\text{yr}+1 -k}$$

This answers the first part of the question in the OP:

Is there an analytical solution to this problem when f,g are constants?

However, there was a second part, which I tried to resolve slop... quickly... leading to too much poetic license, so poorly tolerated in math circles. So after being called on it, that part is now erased, and I'm trying again. First the second part of this question:

What is a suitable numerical method for solving this problem when f,g are arbitrary functions?

Specifically (see initial comments) the functions Alexandre has in mind are linearly changing probabilities over time as:

$\color{blue}{m_t} = f(\text{yr})= \text{max} \{0, 1−a\times\text{yr}\}$


$\color{blue}{o_t} = g(\text{yr})= \text{max} \{0, 1−b\times \text{yr}\}.$

The problem then becomes apparent in, for example, the expression above for $\small \Pr(\text{Off= }k\,\vert \,\text{mate @ Yr=yr})$: treated as a binomial, it prompts to select any combinations of $k$ successes (offspring), and treats each one with fixed probabilities of success and failure. Unfortunately, the summations over years that follow don't correct for this (my oversight).

So back to the drawing board... The hope is that generalizing $\text{Eq. } *$ will do the trick for the rest of the post. This equation is the multiplication of two probabilities. Starting with the least challenging term:

$$\Pr(\text{mate @ Yr=yr})= m_{\small t=\text{yr}}\,\prod_{t=0}^\text{yr}(1-m_t)$$

This is probably clear, although technically, it could be labeled as a geometric distribution with varying probability values.

The challenge, then, is in $\Pr(\text{Off= }k\,\vert \,\text{mate @ Yr=yr})$. This turns out to be a Poisson binomial distribution, and adapting the notation to our case would result in a beautiful expression:

$$\Pr(\text{Off= }k\text{ by time } x\,\vert \,\text{mate @ Yr=yr})=\sum_{A\in F_k}\,\,\prod_{i\in A}\,o_{t\in i}\,\prod_{j\in A^c}(1-o_{t\in j})$$

As in the Wikipedia link, $F_k$ is the set of all subsets of $k$ integers selected from $\{\text{yr}, \text{yr}+1, \cdots, x\}.$

And just when I was searching for a link to The Scream imagining the process of actually inserting this thing into the other two equations, I realized that the actual challenge would be the numerical calculations, at which point @Zen came to save the day, together with @wolfies.

So just for fun, the final equation formulating the probability of a certain number of offspring ($k$) by a given age ($x$, capped at 10 years in the code formulations below), and having found a mate at year $\text{yr}$, would now look something along the lines of:

$$\Pr(\text{Off by time } x\geq k)= \sum_{\text{Yr=0}}^{x-k+1}\,\, \sum_{\text{K}=k}^{x-\text{yr}+1} \,\,\small \left( m_{\small t=\text{yr}}\,\prod_{t=0}^\text{yr}(1-m_t)\right) \left(\sum_{A\in F_k}\,\,\prod_{i\in A}\,o_{t\in i}\,\prod_{j\in A^c}(1-o_{t\in j})\right)$$

Finally, what would this look like in R:

Yr = 1:11 # Years zero to 10. x corresponds to 10, which is Yr[11].

a = .055  
# Arbritarily chosen slope of the function for p(mating at time = yr)   
(m_t = ifelse((1 - a * Yr) > 0, 1 - a * Yr, 0)) 
# Ifelse to reject potential negative prob. values

# [1] 0.945 0.890 0.835 0.780 0.725 0.670 0.615 0.560 0.505 0.450 0.395

b = .09   
# Doing the same for the slope to calculate P(offspring | mate)
(o_t = ifelse((1 - b * Yr) > 0, 1 - b * Yr, 0)) # Same trick to avoid neg values

# [1] 0.91 0.82 0.73 0.64 0.55 0.46 0.37 0.28 0.19 0.10 0.01

Prob_Off_k_and_mate_yr = function(k=1:11, yr= 0:10){ 
  # Probability Offspring = k AND mating at Yr = yr
  #k needs to be between yr and x
  if(k > (length(Yr) - yr + 1)){stop('Number of spring selected is impossible')}else{
    #Probability to mate at year Yr = yr:
    P_mate_at_yr = ifelse(yr==0, m_t[1], prod(1 - (m_t)[1:yr]) * m_t[yr + 1])
    #Probability of Offspring = k having mated at Yr = yr
    S = seq(yr, length(Yr) - 1) 
    # All the years remaning to choose from, including the mating year.
    A = combn(S, k)             
    # All possible combinations of k times from S
    P_off_k_having_mated_yr = 0 # Starting an empty vector
    for (i in 1:ncol(A)) {    
      # For all subsets of k elements from the years "available"
      P_off_k_having_mated_yr = P_off_k_having_mated_yr +
        prod(o_t[A[,i] + 1], 1 - o_t[setdiff(S, A[,i]) + 1]) 
      # Poisson binomial
    Prob_Off_k_and_mate_yr = P_mate_at_yr * P_off_k_having_mated_yr 

# Trying the function for Offspring = 6 and mating at year 1:
k = 6
yr= 1

# [1] 0.005674715

# What about the probability of Offspring = 6 
# regardless of the mating year (summation over years):
Prob_Off_k = 0
for(i in 0:(length(Yr) - k)){
  Prob_Off_k =  Prob_Off_k + Prob_Off_k_and_mate_yr(k, i)

# [1] 0.2238927

# Finally, the actual question in the OP: AT LEAST 3 Offspring (for example):
Prob_at_least_k = 0 # Starting empty vector
for(i in 0:(length(Yr)- k)){     
  # Loop over mating year, which can't go beyond len(Yr) - k
  Prob_Off_k = 0                
  # Probability of k and any max allowable k depending on the year of mating (i)
  for(j in k:(length(Yr) - i)){ # Index for k's
    Prob_Off_k =  Prob_Off_k + Prob_Off_k_and_mate_yr(j, i)
  Prob_at_least_k = Prob_at_least_k + Prob_Off_k

# [1] 0.9682951

A more elegant way of coding this process, thanks, so many thanks to the wisdom of whuber, in this instance on this post, could be achieved using a convolution with the R function convolve(), which calculates it using a Fast Fourier Transform (FFT). This would be the modified Prob_Off_k_and_mate_yr function:

Prob_Off_k_and_mate_yr_convolution = function(k=1:11, yr= 0:10){ 
  # Probability Offspring = k AND mating at Yr = yr
  #k needs to be between yr and x
  if(k > (length(Yr) - yr + 1)){stop('Number of spring selected is impossible')}else{
    #Probability to mate at year Yr = yr:
    P_mate_at_yr = ifelse(yr==0, m_t[1], prod(1 - (m_t)[1:yr]) * m_t[yr + 1])
    #Probability of Offspring = k having mated at Yr = yr

    z = 1
    for (u in sort(o_t[yr+1:length(o_t)])) z <- convolve(z, c(u, 1 - u), type = "open")
    Prob_Offspring = z *  P_mate_at_yr 
    return(Prob_Offspring[k + 1])

It is shorter in coding lines, more elegant mathematically, and so much faster:

microbenchmark(Prob_Off_k_and_mate_yr_convolution(k,yr), Prob_Off_k_and_mate_yr(k,yr))
Unit: microseconds
                                      expr      min        lq      mean    median       uq      max neval
 Prob_Off_k_and_mate_yr_convolution(k, yr)  281.220  288.7165  298.6452  294.5675  301.333  376.300   100
             Prob_Off_k_and_mate_yr(k, yr) 3959.012 4046.9615 4236.5416 4111.1405 4187.023 6195.602   100
  • $\begingroup$ Thank you for the fantastic answer. The only thing missing is the case where $f,g$ are not constants. Can you shed some light on that so that I can accept your answer? $\endgroup$
    – Alexandre
    Sep 28 '16 at 19:48
  • $\begingroup$ @Alexandre Tough... So I can mull over it and feel inadequate :-) ... in plain English the second part of the question makes reference to a case where we still have a Bernoulli experiment in finding a mate with a given probability, and having offspring (provided a mate exists), but with the parameter $p$ in $X \sim \text{Bern}(p)$ changing from year to year? $\endgroup$ Sep 28 '16 at 19:56
  • $\begingroup$ Would the inter-annual change follow any pattern? $\endgroup$ Sep 28 '16 at 19:56
  • $\begingroup$ Yes, your first comment is correct. In my case, it will be a linear function that decreases towards 0, but for the purpose of helping others in the future, I think it would be more interesting to allow arbitrary functions. Sorry for giving you so much work! Loved your answer though. I can't find a fault in it, and I can imagine it being useful to others in the future. $\endgroup$
    – Alexandre
    Sep 28 '16 at 20:12
  • $\begingroup$ @Alexandre It would be really cool to just plug in a $f(\text{time})$ into $o$ and $m$ if available. Otherwise, the summation would become an indicial notation nightmare... $\endgroup$ Sep 28 '16 at 20:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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