Re: [程式] R: 改寫程式 不使用for loop

[celestialgod (天)]

Andrew大給了apply的做法

但apply的本質還是迴圈

我覺得這問題如果是作業的話

最好的解答就是把var的公式展開，直接做計算

如果不是作業，為了開發的話，最好方式就是直接用var就好

var(t(X)) 即可

直接計算的話，你可以先把公式展開:

   var(x, y) = sum_i (X_i - X_bar) (Y_i - Y_bar) / (ncol(X)-1)

where X_bar = the mean of X, Y_bar = the mean of Y

所以根據這個公式你可以擴展到向量版本

var(X, Y)第i列第j行的元素可以寫作下面這樣

    var(X, Y)_ij = sum_k (X_ki - X_bar_i) (Y_kj - Y_bar_j) / (ncol(X)-1)

                 = (X_i - X_bar_i) (Y_j - Y_bar_j)^T / (ncol(X)-1)

where X and Y are P x N matrices, ^T is transpose,

      X_bar_i = the mean of i-th column of X,

      Y_bar_j = the mean of j-th column of Y,

      X_ij = the data vector of i-th row and j-th column of X,

      Y_ij = the data vector of i-th row and j-th column of Y,

       X_i is the data vector of i-th column of X,

       Y_j is the data vector of j-th column of Y.

如果你稍微再簡化一下，你就會得到

     var(X, Y) = (X - X_bar * I_{p}) (Y - Y_bar * I_{p})^T / (ncol(X)-1)

where X_bar is the data vector of means of the rows of X,

      Y_bar is the data vector of means of the rows of Y,

      I_{p} is p x p identity matrix.

最後用程式表示的話

  temp <- sweep(X, 1, rowMeans(X), `-`)
  alpha <- temp %*% t(temp) / (ncol(X)-1)

結果: https://ideone.com/hgruOe

# Unit: milliseconds
#          expr       min        lq      mean    median        uq       max
#      for_loop  791.3887  814.2140  876.8024  872.5140  907.3111 1071.6495
#    apply_func 1486.9081 1542.1674 1607.6525 1592.6940 1664.5603 1796.9124
#   diag_method  802.1050  825.0755  938.7255  888.7893  973.0518 1453.5382
#  formula_calc  105.0659  110.9181  133.3554  118.8704  140.6495  249.5756
#           var  122.8382  124.8452  125.7241  125.1085  126.3662  130.1424

你就會發現用formula直接計算還會比用var少掉不少時間

因為略過了不少檢查(攤手

不過這個大概就是極限了，主要計算瓶頸在matrix計算

這部分由BLAS負責，我使用的BLAS已經是公認最快的Intel MKL了(攤手

你就算用Rcpp也沒辦法增加多少速度

by 剛好休假有空回文的前R板板主XDDDDDD

※ 引述《smellyhead (smellyheadN       )》之銘言：
: [R code]
: X<-matrix(rnorm(300),3,100)
: alpha <- matrix(c(3,2,1,2,2,1,1,1,1),3,3)
: for(i in 1:3){
:   for(j in i:3){
:     alpha[i,j] <- alpha[j,i] <- var(X[i,],X[j,])
: }
: }
: 請問一下
: 上面程式 要如何改寫成不用for loop
: 謝謝!!

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 114.24.106.91 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/Statistics/M.1617770863.A.74C.html
推 andrew43: 推向量化，但對數學的基本功要求高不少。
  但是統計跟寫code本來就是需要大量數學的XDDDD

  不要跟我說數學不好 可以把這兩個都學好(茶

  所以有些現實還是要早一點讓學習者認清

  數學不好，程式跟統計都很難到達top的水準

  雖說不是每個人都需要到top的水準

  但能有機會進步，誰不願意呢XDD

推 joshddd: 推 不過有時候就是懶吧
懶惰才更該把程式寫好（茶

技術債才不會越積越多，後面才不會改一個東西動輒得咎（攤手

再來就是這是寫code時間的花費與計算時間之間的取捨

你花多一點時間寫好的架構跟速度快的程式

跟花少少時間寫能動速度不快的程式

雖說都能得到一樣的結果

但是這套程式 你日後是不是還要再使用、會不會要修改它會有很大的不同

如果今日懶惰，之後又要新增很多功能的時候

後者保證讓你痛不欲生（笑
※ 編輯: celestialgod (101.10.76.194 臺灣), 04/08/2021 09:35:20
推 joshddd: 謝謝大大的建議 我會改進的
→ DIDIMIN: 懶惰才是寫程式的原動力，只想寫出可以一鍵完成的程式