使用应用函数重写循环

Question

我有以下 3 个我想让它更快的函数，我认为应用函数是最好的方法，但我从未使用过应用函数，所以我不知道该怎么做。任何类型的提示、想法和代码片段将不胜感激。

n、T、dt 是全局参数，par 是参数向量。

函数 1:是创建 m+1,n 矩阵的函数，该矩阵包含具有指数分布跳跃大小的泊松分布跳跃。我的麻烦是因为我有 3 个循环，我不确定如何将 if 语句合并到内部循环中。我也不知道是否可以只在循环的外层使用应用函数。

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 to decide when the jumps will happen
  Y <-replicate(n,rexp(100,1/par[6])) #matrix with jump sizes
  for (l in 1:n){ 
    NT <- rpois(1,par[5]*T) #number of jumps
    k=0
    for (j in seq(t,T,dt)){
      k=k+1
      if (NT>0){
        temp=0
        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]
        }
        jump[k,l]=temp
      }else jump[k,l]=0
    }
  }
  return(jump)
}

函数 2:根据布朗运动和函数 1 的跳跃计算默认强度。这里我的问题是如何使用应用函数，当用于计算的变量是输出矩阵中上面行的值以及如何获得计算中使用的外部矩阵的正确值(BMz_C & J)

lambda <- function(t=0,T=T,par,fit=0){
  lambda <- matrix(0,m+1,n) # matrix to hold intesity path output
  lambda[1,] <- par[4] #initializing start value of the 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)
} 

函数 3:根据函数 2 的强度计算生存概率。这里 a() 和 B() 是返回数值的函数。我的问题是同时使用了 i 和 j 值，因为 i 并不总是一个整数，因此可以用来引用外部矩阵。我之前曾尝试使用 i/dt，但有时它会覆盖矩阵中的一行并跳过下一行，很可能是由于舍入错误。

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
  for (i in seq(t,T,dt)){
    j=j+1
    S[j,] <- a(i,T,par)*exp(B(i,T,par)*l[j,])
  }
  return(S)
}

很抱歉这篇长文章，对任何功能的任何帮助将不胜感激。

编辑:
首先感谢 digEmAll 的精彩回复。

我现在已经研究了矢量化功能 2。首先我尝试过

lambda <- function(t=0,T=T,par,fit=0){
  lambda <- matrix(0,m+1,n) # matrix to 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)
}

但它只会产生第一列。所以我尝试了一个两步应用函数。

lambda <- function(t=0,T=T,par,fit=0){
  lambda <- matrix(0,m+1,n) # matrix to 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)
}

这似乎有效，但仅在第一行，之后的所有行都具有相同的非零值，就好像 lambda[i-1] 未用于 lambda[i] 的计算一样，有没有人知道如何管理那个？

Answer 1

我将逐步向您解释如何矢量化第一个函数( 一种可能的矢量化方式 ，可能不是最适合您的情况)。
对于其他 2 个功能，您可以简单地应用相同的概念，您应该能够做到。

这里的关键概念是:从最内层循环开始向量化。

1) 首先，rpois一次可以生成一个以上的随机值，但您调用它 n 次询问一个随机值。所以，让我们把它从获得这个的循环中取出来:

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 the change
  for (l in 1:n){ 
    NT <- NTs[l] # note the change
    k=0
    for (j in seq(t,T,dt)){
      k=k+1
      if (NT>0){
        temp=0
        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]
        }
        jump[k,l]=temp
      }else jump[k,l]=0
    }
  }
  return(jump)
}

2) 同样，调用 seq(t,T,dt) 没用/效率低下循环中的 n 次，因为它总是会生成相同的序列。所以，让我们把它从循环中取出并存储到一个向量中，得到这个:

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 the change
  for (l in 1:n){ 
    NT <- NTs[l] 
    k=0
    for (j in js){ # note the change
      k=k+1
      if (NT>0){
        temp=0
        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]
        }
        jump[k,l]=temp
      }else jump[k,l]=0
    }
  }
  return(jump)
}

3) 现在，让我们看看最里面的循环:

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]
}

这等于:

u <- as.integer(U[1:NT,l]<=j)
temp <- sum(Y[1:NT,l]*u)

或者，在一行中:

temp <- sum(Y[1:NT,l] * as.integer(U[1:NT,l] <= j))

因此，现在该函数可以写为:

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) 
  for (l in 1:n){ 
    NT <- NTs[l] 
    k=0
    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 the change
      }else jump[k,l]=0
    }
  }
  return(jump)
}

4) 再一次，让我们看看当前最里面的循环:

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 the change
  }else jump[k,l]=0
}

如您所见， NT不依赖于这个循环的迭代，所以内部 if可以移到外面，如下:

if (NT>0){
  for (j in js){ 
    k=k+1
    jump[k,l] <- sum(Y[1:NT,l]*as.integer(U[1:NT,l]<=j)) # note the change
  }
}else{
  for (j in js){ 
    k=k+1
    jump[k,l]=0
  }
}

这似乎比以前更糟，是的，但现在这两个条件可以变成单行的(注意使用 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
}

获得如下跳转函数:

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) 
  for (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) 最后我们可以摆脱最后一个循环，再次使用 sapply¹函数，获得最终 jump功能 :

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 功能很容易使用。对于在 X 中传递的列表或向量的每个元素参数，它应用在 FUN 中传递的函数参数，例如:

vect <- 1:3
sapply(X=vect,FUN=function(el){el+10}
# [1] 11 12 13

因为默认情况下 simplify参数为真，结果被强制为最简单的可能对象。因此，例如在前面的例子中结果变成了一个向量，而在下面的例子中结果变成了一个矩阵(因为对于每个元素我们返回一个相同大小的向量):

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

基准:

以下基准测试只是让您了解速度增益，但实际性能可能会因您的输入参数而异。
可以想象， jump_old对应于你原来的函数1，而 jump_new是最终的矢量化版本。

# let's use some 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 the 2 functions are the same:
isTRUE(all.equal(old,new))
# [1] TRUE