Faster alternative to for loops in R?

user_10161382 0
Viewed 97

I am trying to calculate a transition probability (without Markov assumption), which require calculating this nested integration, nested integral

Note that the integrals can be replace by summations in my case. Here is a toy example code I am using to calculate this,

# simulate some data
set.seed(99)
data<-data.frame(time=seq(0,7,0.1),
             S_D=seq(1,0.95,length.out = 71),
             lam12=sample(c(0,0.1,0.12,0.15,0.17),size = 71,replace = TRUE),
             lam23=sample(c(0,0.05,0.1,0.08,0.12),size = 71,replace = TRUE),
             lam24=sample(c(0,0.02,0.05,0.06,0.08),size = 71,replace = TRUE))

prob_123<-c() # initializing a NULL vector
end<-nrow(data)

for (j in 2: end)
{
  # j indicates u in the expresstion
  # k indicates v in the expression
  prob_123k<-0
  for (k in (j+1):end)
  {
    if (k==(j+1)){  
      prob_123k<-prob_123k+data$S_D[j-1]*data$lam12[j]*data$lam12[k-j]
    }
    if (k>(j+1)){ 
  
      prob_123k<-prob_123k+data$S_D[j-1]*data$lam12[j]*prod(1-(data$lam12[1:(k-j-1)]+data$lam24[1:(k-j-1)]))*(data$lam12[k-j])
   }
  }
  prob_123[j-1]<-prob_123k
}

sum(prob_123)  # result = 5.631623

In the code, S_D corresponds to expression exp{-(\Lambda_12(u)+\Lambda_13(u)+\Lambda_14(u))} and prod(1-(...)) corresponds to expression exp{-(\Lambda_23(v-u)+\Lambda_24(v-u))}. My original dataset is much bigger than this and it takes a long time calculate the nested for loops. Can anyone please suggest any faster alternatives? Thank you so much.

r loops
1 Answers

You could use a combination of cumprod and sum to get the loops faster.

Encapsulating your for-loop into a function f:

f <- function(data) {
  end <- nrow(data)
  prob_123 <- vector("numeric", end) # initializing a NULL vector
  for (j in 2:end) {
    # j indicates u in the expresstion
    # k indicates v in the expression
    prob_123k <- 0
    for (k in (j + 1):end) {
      
      if (k == (j + 1)) {  
        prob_123k <- prob_123k + data$S_D[j - 1] * data$lam12[j] * data$lam12[k - j]
      }
      if (k > (j + 1)) { 
        prob_123k <- prob_123k + data$S_D[j - 1] * data$lam12[j] * data$lam12[k - j] *
          prod(1 - (data$lam12[1:(k - j - 1)] + data$lam24[1:(k - j - 1)]))
      }
    }
    prob_123[j - 1] <- prob_123k
  }
  
  sum(prob_123)    
}

f(data)
5.631623

We note that there's an unnecessary re-computation of the running product - it can be computed once and then indexed appropriately. Rewriting the function as follows:

ff <- function(data) {
  end <- nrow(data)
  prob_123 <- vector("numeric", end) # initializing a NULL vector
  p <- cumprod(1 - (data$lam12 + data$lam24))
  for (j in 2:end) {
    # j indicates u in the expresstion
    # k indicates v in the expression
    if (j + 1 <= end) {
      prob_123k <- data$S_D[j - 1] * data$lam12[j] * data$lam12[(j + 1):end - j] * 
        c(p[1], p[1:(end - j - 1)])
    }
    prob_123[j - 1] <- sum(prob_123k) + # Equivalent to the k > j + 1 part
      data$S_D[j - 1] * data$lam12[j] * data$lam12[j + 1 - j] # Equivalent to k = j + 1 part
  }
  
  sum(prob_123)    
}

identical(f(data), ff(data))
TRUE

We can then use the microbenchmark package to see if there's an improvement:

library(microbenchmark)
microbenchmark(f(data), ff(data))
Unit: microseconds
     expr      min       lq     mean   median        uq       max neval
  f(data) 8853.501 9020.902 10118.14 9410.351 10416.201 14392.401   100
 ff(data)  344.701  356.401   373.86  367.301   384.701   466.401   100

The function ff is about 30x faster on average.

I'm sure code jockeys can optimise this even more, by getting rid of j for-loop, perhaps?

user_24742260
Related