【发布时间】：2014-07-21 18:30:26
【问题描述】：
我有以下 3 个想要更快的函数，我认为应用函数是最好的方法，但我从未使用过应用函数，所以我不知道该怎么做。任何类型的提示、想法和代码 sn-ps 将不胜感激。
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]，有没有人有知道如何管理它？
【问题讨论】：
-
旁注：你知道
T等于True吗？您真的应该考虑更改该参数的名称以避免混淆...... -
如果您从未使用过
apply系列函数，您可以先看看this question -
您应该发布一个调用函数的示例（带有有效参数和全局参数的示例）。通过这种方式，我们将有一个可重现的示例来处理......
-
@Marcus：同样，如果您希望我们能够帮助您，您需要发布输入数据的示例（或者可能是较小的子集）。我的意思是：
m,n,dt,T,par,BMz_C等等 ... -
@Marcus：顺便说一句，我想补充一句。您使用我的代码获得的速度提升主要是由于在循环外提取了
rpois和seq调用。事实上，sapply并没有比写得好的循环快 太多（“写得好”我的意思是一个你不附加到矩阵而不是预分配的循环，你不t 多次重新计算可以存储和重复使用的东西等）。另外，你的第二个乐趣。似乎需要上一次迭代的结果，这很难使用 sapply 进行矢量化，因此，如果性能可以接受，我会坚持使用你的循环；）
标签： r for-loop recursion apply performance