0
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$
    – jbowman
    Jul 21, 2017 at 13:26

1 Answer 1

1
$\begingroup$

It looks like recursively calling Prob_S is killing you. This looks like a classic case where dynamic programming may be used. The idea is simple- don't do the same calculation twice. Keep a table where results of every call to Prob_S is stored and when calling Prob_S, check that table to see if it's a calculation you've done before.

A classic example is writing a recursive program to find the nth Fibonacci number. If you recurse blindly you do the same computation again and again. However, if you use dynamic programming and store the results of your previous computations, you get a significant speedup. Check Wikipedia, the Fibonacci number calculation example may be of use.

$\endgroup$
4
  • $\begingroup$ I agree with your diagnosis, but it doesn't look to me as though there are any repeats of the calls to Prob_S (in the sense of the same parameters being passed), unlike what happens in a Fibonacci recursion. If this is correct, DP / memoization won't help, or so it seems to me. $\endgroup$
    – jbowman
    Jul 21, 2017 at 14:00
  • $\begingroup$ Dear Yair Daon, you were right. I solved it. Thanks Yair Daon and jbowman for your help. Kind regards. $\endgroup$ Sep 3, 2017 at 4:56
  • $\begingroup$ Glad you solved it. So was it dynamic programming or something else, as @jbowman suggested? $\endgroup$
    – Yair Daon
    Sep 3, 2017 at 6:41
  • $\begingroup$ I used dynamic programming and the efficiency of the algorithm improved. $\endgroup$ Sep 3, 2017 at 11:43

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.