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}\}$
and
$\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
return(Prob_Off_k_and_mate_yr)
}
}
# Trying the function for Offspring = 6 and mating at year 1:
k = 6
yr= 1
Prob_Off_k_and_mate_yr(k,yr)
# [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)
}
Prob_Off_k
# [1] 0.2238927
# Finally, the actual question in the OP: AT LEAST 3 Offspring (for example):
k=3
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
}
Prob_at_least_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
