# How to compute more efficiently in R the probability distribution of the sum of non-independent discrete random variables

I hope you are well.
Let $\{s_0,\,s_1,\ldots,\,s_T\}$ be a sequence of discrete random variables and denote $S_t=s_0+s_1+\cdots+s_t$, with $S_0=0$.
For all $t\in\{1,\ldots,\,T\}$, suppose that
$s_t|\{S_{t-1},\,p_t\}\sim\text{Binomial}(M-S_{t-1},\,p_t)$, with $M$ is a fixed positive integer,
$\text{logit}(p_t)=\beta_0+\beta_1\cdot S_{t-1}$, and $\beta_0\in\mathbb{R}$ and $\beta_1\in\mathbb{R}$ are known and fixed.
Conditionally on $M$, $\beta_0$, $\beta_1$, $t$, and $m$, with $m\in\{0,\,1,\ldots,\,M\}$, the following R-code is intended to compute $\mathbb{P}(S_t=m)$.

library(boot)

TMax <- 20        # In this R-code, I am using TMax instead of using T.
M <- 100
beta0 <- 1
beta1 <- 0.5
Prob_S <- function(m, r){        # In this R-code, I am using r instead of using t.
if(r == 1){
Aux <- dbinom(x = m, size = M, prob = inv.logit(beta0))
}
if(r %in% 2:TMax){
Aux <- 0
for(u in 0:m){
Aux <- Aux + dbinom(x = m - u, size = M - u,
prob = inv.logit(beta0 + beta1 * u)) * Prob_S(u, r - 1)
}
}
Aux
}


This R-code builds on the recursive formula:
$\displaystyle\mathbb{P}(S_1=k_1)={M\choose k_1}\cdot p_1^{k_1}\cdot(1-p_1)^{M-k_1}$ and

$\displaystyle\mathbb{P}(S_t=k_t)=\sum_{k_{t-1}=0}^{k_t}\mathbb{P}(s_t=k_t-k_{t-1}|\{S_{t-1}=k_{t-1}\})\cdot\mathbb{P}(S_{t-1}=k_{t-1})$,
for all $t\in\{2,\ldots,\,T\}$. However, I realized that this R-code is inefficient. For example, the command

Prob_S(m = 15, r = 10)


takes several hours to compute $\mathbb{P}(S_{10}=15)$.
Question: How can I compute $\mathbb{P}(S_t=m)$ more efficiently in R?
Thanks a lot for your help and suggestions.

• Just out of curiosity, what sort of computer are you running on? It takes 24 seconds (still worthy of investigation of course) on my 9 year old 8 core 6GB Windows desktop, R version 3.4.x. Jul 21, 2017 at 13:26