用apply函数重写循环

Question

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

n,T,dt是全局参数,par是参数的矢量.

功能1:是一个创建m + 1,n矩阵的函数,该矩阵包含具有指数分布的跳跃大小的泊松分布跳跃.我的麻烦是因为我有3个循环,我不知道如何在内循环中加入if语句.另外我不知道是否完全可以在循环的外层使用apply函数.

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参数为true,结果将被强制转换为最简单的对象.因此,例如在前一种情况下,结果变为向量,而在下面的示例中,结果变为矩阵(因为对于每个元素,我们返回相同大小的向量):

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