计算总和嵌套环中的纵向数据不平衡

Question

我需要按下嵌套环中的每个i上的每个j上的总和，然后输出每个i作为列表的总和。问题在于，对于大量观察，代码的速度非常慢。有什么方法可以避免循环以使代码运行速度更快？谢谢。

#### generate data
 set.seed(234)
 N=3
 v<-sample(2:6,N,replace=TRUE)
 id<-c(rep(1:N,v))
 n<-length(id)
 x<-as.matrix(cbind(rnorm(n,0,1),rnorm(n,0,1),rnorm(n,0,1)))
 x1<-cbind(id,x)
 e<-runif(3)
 > v
 [1] 5 5 2
 id
  [1] 1 1 1 1 1 2 2 2 2 2 3 3
> x
            [,1]       [,2]        [,3]
 [1,]  0.7590390 -0.8716028 -0.30554099
 [2,]  0.3713058  1.1876234  0.86956546
 [3,]  0.5758514 -0.6672287 -1.06121591
 [4,] -0.5703207  0.5383396 -0.09635967
 [5,]  0.1198567  0.4905632  0.47460932
 [6,]  0.2095484 -1.0216529 -0.02671707
 [7,] -0.1481357 -0.3726091  1.10167492
 [8,]  0.6433900  1.3251178 -0.26842418
 [9,]  1.1348350 -0.7313432  0.01035965
[10,]  0.1995994  0.7625386  0.25897152
[11,]  0.2987197  0.3275333 -0.39459737
[12,] -0.3191671 -1.1440187 -0.48873668
> e
[1] 0.3800745 0.5497359 0.3893235

 ### compute sum
  sumterm_<-list()
  count=1
 for (i in 1:N){
   idd=x1[,1]==i
   xi=x[idd,]
  sumterm=matrix(rep(0,N*N),nrow=3,ncol=3)
  for (j in 1:v[i]){
    xij=xi[j,]
     sumterm=sumterm+as.matrix(xij-e)%*%(xij-e)
     count=count+1
  }
   sumterm_[[i]]<-sumterm
  }
sumterm_
[[1]]
           [,1]       [,2]       [,3]
[1,]  1.1529838 -0.7562553 -0.1121242
[2,] -0.7562553  3.9117383  3.0597216
[3,] -0.1121242  3.0597216  3.0606953
[[2]]
             [,1]        [,2]        [,3]
 [1,]  0.97965490 -0.04598867 -0.74102232
 [2,] -0.04598867  5.60764839 -0.05553464
 [3,] -0.74102232 -0.05553464  1.27377151
[[3]]
          [,1]     [,2]      [,3]
[1,] 0.4955573 1.202421 0.6777518
[2,] 1.2024208 2.918179 1.6614076
[3,] 0.6777518 1.661408 1.3855215

Answer 1

可以采取一些步骤来改进代码：

• 一个GO

  为输出对象分配所有空间

  sumterm_ <- lapply(1:N,function(x){matrix(0,3,3)})

• 一次计算X-E一次，而不是重复相同的二元

  xbar <- x-rep(e, each=n)

• 使用drop=FALSE避免将矩阵转换为向量并返回

  xbar[i,] %*% xbar[i,,drop=FALSE]

• 直接写入输出对象

  sumterm_[[id[i]]] <- sumterm_[[id[i]]] + xbar[i,] %*% xbar[i,,drop=FALSE]

因此，完整的代码看起来像：

  #List of zero matrices
  sumterm_ <- lapply(1:N,function(x){matrix(0,3,3)})
  #Calculate x-e
  xbar <- x-rep(e, each=n)
  #sum by id 
  for (i in 1:n){
    sumterm_[[id[i]]] <- sumterm_[[id[i]]] + xbar[i,] %*% xbar[i,,drop=FALSE]
  }

另一种方法可能是使用应用功能重写(尽管这些功能在其中实现了，而不是消除它们(。

#calculate cross product for each row
cps <- apply(x-rep(e, each=n), 1, tcrossprod)
#aggregate columns by id variable, and convert to matrix
sumterm2_ <- tapply(seq_along(id), id, 
                    function(i){matrix(rowSums(cps[, i, drop=FALSE]), 3, 3)})

比较不同方法之间的速度取决于问题范围的方向 - 这就是为什么方法之间没有时间比较的原因。