title: Rewriting loops with apply functions
body: <p>I have the 3 following functions which I would like to make faster, I assume apply functions are the best way to go, but I have never used apply functions, so I have no idea what to do. Any type of hints, ideas and code snippets will be much appreciated. </p>
  
  <p>n, T, dt are global parameters and par is a vector of parameters.</p>
  
  <p>Function 1: is a function to create an m+1,n matrix containing poisson distributed jumps with exponentially distributed jump sizes. My troubles here is because I have 3 loops and I am not sure how to incorporate the if statement in the inner loop. Also I have no idea if it is at all possible to use apply functions on the outer layers of the loops only.</p>
  
  <pre><code>jump &lt;- function(t=0,T=T,par){
    jump &lt;- matrix(0,T/dt+1,n) # initializing output matrix
    U &lt;- replicate(n,runif(100,t,T)) #matrix used to decide when the jumps will happen
    Y &lt;-replicate(n,rexp(100,1/par[6])) #matrix with jump sizes
    for (l in 1:n){ 
      NT &lt;- rpois(1,par[5]*T) #number of jumps
      k=0
      for (j in seq(t,T,dt)){
        k=k+1
        if (NT&gt;0){
          temp=0
          for (i in 1:NT){
            u &lt;- vector("numeric",NT)
            if (U[i,l]&gt;j){ u[i]=0
            }else u[i]=1
            temp=temp+Y[i,l]*u[i]
          }
          jump[k,l]=temp
        }else jump[k,l]=0
      }
    }
    return(jump)
  }
  </code></pre>
  
  <p>Function 2: calculates a default intensity, based on Brownian motions and the jumps from function 1. Here my trouble is how to use apply functions when the variable used for the calculation is the values from the row above in the output matrix AND how to get the right values from the external matrices which are used in the calculations (BMz_C &amp; J)</p>
  
  <pre><code>lambda &lt;- function(t=0,T=T,par,fit=0){
    lambda &lt;- matrix(0,m+1,n) # matrix to hold intesity path output
    lambda[1,] &lt;- par[4] #initializing start value of the intensity path.
    J &lt;- jump(t,T,par) #matrix containing jumps
    for(i in 2:(m+1)){
      dlambda &lt;- par[1]*(par[2]-max(lambda[i-1,],0))*dt+par[3]*sqrt(max(lambda[i-   1,],0))*BMz_C[i,]+(J[i,]-J[i-1,])
      lambda[i,] &lt;- lambda[i-1,]+dlambda
    } 
    return(lambda)
  } 
  </code></pre>
  
  <p>Function 3: calculates a survival probability based on the intensity from function 2. Here a() and B() are functions that return numerical values. My problem here is that the both value i and j are used because i is not always an integer which thus can to be used to reference the external matrix. I have earlier tried to use i/dt, but sometimes it would overwrite one line and skip the next lines in the matrix, most likely due to rounding errors.</p>
  
  <pre><code>S &lt;- function(t=0,T=T,par,plot=0, fit=0){
    S &lt;- matrix(0,(T-t)/dt+1,n)
    if (fit &gt; 0) S.fit &lt;- matrix(0,1,length(mat)) else S.fit &lt;- 0
    l=lambda(t,T,par,fit)
    j=0
    for (i in seq(t,T,dt)){
      j=j+1
      S[j,] &lt;- a(i,T,par)*exp(B(i,T,par)*l[j,])
    }
    return(S)
  }
  </code></pre>
  
  <p>Sorry for the long post, any help for any of the functions will be much appreciated.</p>
  
  <p>EDIT:
  First of all thanks to digEmAll for the great reply.</p>
  
  <p>I have now worked on vectorising function 2. First I tried </p>
  
  <pre><code>lambda &lt;- function(t=0,T=T,par,fit=0){
    lambda &lt;- matrix(0,m+1,n) # matrix to hold intesity path input
    J &lt;- jump(t,T,par,fit)
      lambda[1,] &lt;- par[4] 
      lambda[2:(m+1),] &lt;- sapply(2:(m+1), function(i){
        lambda[i-1,]+par[1]*(par[2]-max(lambda[i-1,],0))*dt+par[3]*sqrt(max(lambda[i-1,],0))*BMz_C[i,]+(J[i,]-J[i-1,])
      })
    return(lambda)
  }
  </code></pre>
  
  <p>but it would only produce the first column. So I tried a two step apply function.</p>
  
  <pre><code>lambda &lt;- function(t=0,T=T,par,fit=0){
    lambda &lt;- matrix(0,m+1,n) # matrix to hold intesity path input
    J &lt;- jump(t,T,par,fit)
    lambda[1,] &lt;- par[4] 
    lambda[2:(m+1),] &lt;- sapply(1:n, function(l){
      sapply(2:(m+1), function(i){
        lambda[i-1,l]+par[1]*(par[2]-max(lambda[i-1,l],0))*dt+par[3]*sqrt(max(lambda[i-1,l],0))*BMz_C[i,l]+(J[i,l]-J[i-1,l])
      })
    })
    return(lambda)
  }
  </code></pre>
  
  <p>This seems to work, but only on the first row, all rows after that have an identical non-zero value, as if lambda[i-1] is not used in the calculation of lambda[i], does anyone have an idea how to manage that?</p>