r - Rewriting loops with apply functions -
i have 3 following functions make faster, assume apply functions best way go, have never used apply functions, have no idea do. type of hints, ideas , code snippets appreciated.
n, t, dt global parameters , par vector of parameters.
function 1: function create m+1,n matrix containing poisson distributed jumps exponentially distributed jump sizes. troubles here because have 3 loops , not sure how incorporate if statement in inner loop. have no idea if @ possible use apply functions on outer layers of loops only.
jump <- function(t=0,t=t,par){ jump <- matrix(0,t/dt+1,n) # initializing output matrix u <- replicate(n,runif(100,t,t)) #matrix used decide when jumps happen y <-replicate(n,rexp(100,1/par[6])) #matrix jump sizes (l in 1:n){ nt <- rpois(1,par[5]*t) #number of jumps k=0 (j in seq(t,t,dt)){ k=k+1 if (nt>0){ temp=0 (i in 1:nt){ u <- vector("numeric",nt) if (u[i,l]>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) } function 2: calculates default intensity, based on brownian motions , jumps function 1. here trouble how use apply functions when variable used calculation values row above in output matrix , how right values external matrices used in calculations (bmz_c & j)
lambda <- function(t=0,t=t,par,fit=0){ lambda <- matrix(0,m+1,n) # matrix hold intesity path output lambda[1,] <- par[4] #initializing start value of intensity path. j <- jump(t,t,par) #matrix containing jumps for(i in 2:(m+1)){ dlambda <- 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,] <- lambda[i-1,]+dlambda } return(lambda) } function 3: calculates survival probability based on intensity function 2. here a() , b() functions return numerical values. problem here both value , j used because not integer can used reference external matrix. have earlier tried use i/dt, overwrite 1 line , skip next lines in matrix, due rounding errors.
s <- function(t=0,t=t,par,plot=0, fit=0){ s <- matrix(0,(t-t)/dt+1,n) if (fit > 0) s.fit <- matrix(0,1,length(mat)) else s.fit <- 0 l=lambda(t,t,par,fit) j=0 (i in seq(t,t,dt)){ j=j+1 s[j,] <- a(i,t,par)*exp(b(i,t,par)*l[j,]) } return(s) } sorry long post, any of functions appreciated.
edit: first of digemall great reply.
i have worked on vectorising function 2. first tried
lambda <- function(t=0,t=t,par,fit=0){ lambda <- matrix(0,m+1,n) # matrix hold intesity path input j <- jump(t,t,par,fit) lambda[1,] <- par[4] lambda[2:(m+1),] <- 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) } but produce first column. tried 2 step apply function.
lambda <- function(t=0,t=t,par,fit=0){ lambda <- matrix(0,m+1,n) # matrix hold intesity path input j <- jump(t,t,par,fit) lambda[1,] <- par[4] lambda[2:(m+1),] <- 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) } this seems work, on first row, rows after have identical non-zero value, if lambda[i-1] not used in calculation of lambda[i], have idea how manage that?
i'm going explain you, setp-by-step, how vectorize first function (one possible way of vectorization, maybe not best 1 case).
others 2 functions, can apply same concepts , should able it.
here, key concept is: start vectorize innermost loop.
1) first of all, rpois can generate more 1 random value @ time calling n-times asking 1 random value. so, let's take out of loop obtaining this:
jump <- function(t=0,t=t,par){ jump <- matrix(0,t/dt+1,n) u <- replicate(n,runif(100,t,t)) y <-replicate(n,rexp(100,1/par[6])) nts <- rpois(n,par[5]*t) # note change (l in 1:n){ nt <- nts[l] # note change k=0 (j in seq(t,t,dt)){ k=k+1 if (nt>0){ temp=0 (i in 1:nt){ u <- vector("numeric",nt) if (u[i,l]>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) } 2) similarly, useless/inefficient call seq(t,t,dt) n-times in loop since generate same sequence. so, let's take out of loop , store vector, obtainig this:
jump <- function(t=0,t=t,par){ jump <- matrix(0,t/dt+1,n) u <- replicate(n,runif(100,t,t)) y <-replicate(n,rexp(100,1/par[6])) nts <- rpois(n,par[5]*t) js <- seq(t,t,dt) # note change (l in 1:n){ nt <- nts[l] k=0 (j in js){ # note change k=k+1 if (nt>0){ temp=0 (i in 1:nt){ u <- vector("numeric",nt) if (u[i,l]>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) } 3) now, let's have @ innermost loop:
for (i in 1:nt){ u <- vector("numeric",nt) if (u[i,l]>j){ u[i]=0 }else u[i]=1 temp=temp+y[i,l]*u[i] } this equal :
u <- as.integer(u[1:nt,l]<=j) temp <- sum(y[1:nt,l]*u) or, in one-line:
temp <- sum(y[1:nt,l] * as.integer(u[1:nt,l] <= j)) hence, function can written :
jump <- function(t=0,t=t,par){ jump <- matrix(0,t/dt+1,n) u <- replicate(n,runif(100,t,t)) y <-replicate(n,rexp(100,1/par[6])) nts <- rpois(n,par[5]*t) js <- seq(t,t,dt) (l in 1:n){ nt <- nts[l] k=0 (j in js){ k=k+1 if (nt>0){ jump[k,l] <- sum(y[1:nt,l]*as.integer(u[1:nt,l]<=j)) # note change }else jump[k,l]=0 } } return(jump) } 4) again, let's have @ current innermost loop:
for (j in js){ k=k+1 if (nt>0){ jump[k,l] <- sum(y[1:nt,l]*as.integer(u[1:nt,l]<=j)) # note change }else jump[k,l]=0 } as can notice, nt not depend on iteration of loop, inner if can moved outside, follows:
if (nt>0){ (j in js){ k=k+1 jump[k,l] <- sum(y[1:nt,l]*as.integer(u[1:nt,l]<=j)) # note change } }else{ (j in js){ k=k+1 jump[k,l]=0 } } this seems worse before, yes is, 2 conditions can turned one-liner's (note use of sapply¹):
if (nt>0){ jump[1:length(js),l] <- sapply(js,function(j){ sum(y[1:nt,l]*as.integer(u[1:nt,l]<=j)) }) }else{ jump[1:length(js),l] <- 0 } obtaining following jump function:
jump <- function(t=0,t=t,par){ jump <- matrix(0,t/dt+1,n) u <- replicate(n,runif(100,t,t)) y <-replicate(n,rexp(100,1/par[6])) nts <- rpois(n,par[5]*t) js <- seq(t,t,dt) (l in 1:n){ nt <- nts[l] if (nt>0){ jump[1:length(js),l] <- sapply(js,function(j){ sum(y[1:nt,l]*as.integer(u[1:nt,l]<=j)) }) }else{ jump[1:length(js),l] <- 0 } } return(jump) } 5) can rid of last loop, using again sapply¹ function, obtaining final jump function :
jump <- function(t=0,t=t,par){ u <- replicate(n,runif(100,t,t)) y <-replicate(n,rexp(100,1/par[6])) js <- seq(t,t,dt) nts <- rpois(n,par[5]*t) jump <- sapply(1:n,function(l){ nt <- nts[l] if (nt>0){ sapply(js,function(j){ sum(y[1:nt,l]*as.integer(u[1:nt,l]<=j)) }) }else { rep(0,length(js)) } }) return(jump) } (¹)
sapply function pretty easy use. each element of list or vector passed in x parameter, applies function passed in fun parameter, e.g. :
vect <- 1:3 sapply(x=vect,fun=function(el){el+10} # [1] 11 12 13 since default simplify parameter true, result coerced simplest possible object. so, example in previous case result becomes vector, while in following example result become matrix (since each element return vector of same size) :
vect <- 1:3 sapply(x=vect,fun=function(el){rep(el,5)}) # [,1] [,2] [,3] # [1,] 1 2 3 # [2,] 1 2 3 # [3,] 1 2 3 # [4,] 1 2 3 # [5,] 1 2 3 benchmark :
the following benchmark give idea of speed gain, actual performances may different depending on input parameters.
can imagine, jump_old corresponds original function 1, while jump_new final vectorized version.
# let's use random parameters n = 10 m = 3 t = 13 par = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6) dt <- 3 set.seed(123) system.time(for(i in 1:5000) old <- jump_old(t=t,par=par)) # user system elapsed # 12.39 0.00 12.41 set.seed(123) system.time(for(i in 1:5000) new <- jump_new(t=t,par=par)) # user system elapsed # 4.49 0.00 4.53 # check if last results of 2 functions same: istrue(all.equal(old,new)) # [1] true
Comments
Post a Comment