La probabilidad de apareamiento en cualquier año dado es $\small \Pr(\text{mate})=m$, y la probabilidad de que la descendencia dado a un compañero que se ha encontrado es $\small \Pr( \text{single offspring} \vert \text{mate} ) =o$ y no cambia después de que el compañero se encuentra.
Las probabilidades de contraer $k$ descendencia después de $x$ años depende del año en que la pareja entró en la foto.
Si queremos marcar el comienzo de la "experimento" como el año cero, $\small \text{Yr}=0$, la probabilidad de tener una sola cría al año $0$ es simplemente $\small \Pr(\text{mate } \cap \text{ offspring})= \Pr(\text{offspring}\vert \text{mate})\Pr(\text{mate})=o\times m.$
En el año $1$ (un año después) la probabilidad de que una sola cría, es que simplemente va a ser $\small\Pr(\text{offspring}\vert \text{mate})=o$: la participación de un compañero que ahora está garantizado.
Tan sólo tenemos que tratar cada año después de la aparición de la pareja como un binomio, centrándose en el año en que el mate se "dio cuenta": la probabilidad de que año $1$ es el primer año con una relación de posición se calcula en base a la probabilidad de ausencia de su compañero en el año cero, $1-m$$\small \Pr_{\text{Yr=1}}(\text{mate})=(1-m)\,m$. En el año dos, $\small \Pr_{\text{Yr=2}}(\text{mate})=(1-m)^2\,m$; y, en general, $\small \Pr_{\text{Yr=yr}}(\text{mate})=(1-m)^{\text{yr}}\,m.$
Así tenemos que la probabilidad de tener $k$ crías por año$x$, y que el compañero aparece en el año $\small \text{Yr=yr}$ es:
\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}
Ahora bien, la probabilidad de contraer $k$ crías por año $x$ puede ocurrir en cualquiera de estos escenarios donde la pareja aparece en algunos años o de otro, así:
$$\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{al menos }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}$$
Esto responde a la primera parte de la pregunta del OP:
Existe una solución analítica a este problema si f,g son
constantes?
Sin embargo, hubo una segunda parte, que he tratado de resolver bodrio... rápidamente... que conduce a un exceso de licencia poética, tan mal tolerada en matemáticas círculos. Así que después de ser llamado, esa parte ha sido borrado, y estoy intentando de nuevo. Primero la segunda parte de esta pregunta:
¿Qué es un método numérico para la solución de este problema cuando
f,g son funciones arbitrarias?
Específicamente (ver comentarios iniciales) las funciones Alexandre tiene en mente son linealmente probabilidades que cambian a lo largo del tiempo como:
$\color{blue}{m_t} = f(\text{yr})= \text{max} \{0, 1−a\times\text{yr}\}$
y
$\color{blue}{o_t} = g(\text{yr})= \text{max} \{0, 1−b\times \text{yr}\}.$
El problema entonces se hace evidente, por ejemplo, en la expresión de arriba para $\small \Pr(\text{Off= }k\,\vert \,\text{mate @ Yr=yr})$: tratada como un binomio, se le pedirá que seleccione cualquier combinación de $k$ éxitos (descendencia), y trata a cada uno de los fijos de las probabilidades de éxito y de fracaso. Por desgracia, los sumarios a través de los años que siguen no corregir esto (mi supervisión).
Así que de vuelta a la mesa de dibujo... La esperanza es que la generalización de $\text{Eq. } *$ hará el truco para el resto de la entrada. Esta ecuación es la multiplicación de dos probabilidades. Empezando por el menos difícil plazo:
$$\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{compañero @ 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{año}, \text{año}+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)$$
Por último, ¿qué sería de este aspecto en 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
Una forma más elegante de la codificación de este proceso, gracias, muchas gracias a la sabiduría de whuber, en este caso en este post, podría lograrse mediante una convolución con la función de R convolve(), que se calcula utilizando una transformada Rápida de Fourier (FFT). Esta sería la modificación de la Prob_Off_k_and_mate_yr función de:
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])
}
}
Es más corto en la codificación de líneas, más elegante matemáticamente, y mucho más rápido:
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
