2017-04-24 61 views
1

我需要爲每個嵌套循環中的每個i計算一個總和,如下所示,然後將每個i的總和作爲列表輸出。問題是代碼對於大量的觀察來說非常慢。有什麼辦法可以避免循環,以便代碼運行速度更快?謝謝。在R中計算嵌套循環不平衡縱向數據

#### 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 
+0

你正在成長你的對象。 R中的菜鳥101錯誤預先分配所有內容並填寫結果。 。請參閱[R Inferno](http://www.burns-stat.com/documents/books/the-r-inferno/)獲取有關如何避免不當行爲的更多提示。 –

回答

1

可以採取提高代碼的一些步驟:

  • 分配所有的空間,爲您輸出對象在一個去

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

  • 計算XE一次,而不是重複相同的計算

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

比較不同方法之間的速度取決於問題擴展的方向 - 這就是爲什麼方法之間沒有時間比較的原因。