Statistics 板


LINE

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://webptt.com/m.aspx?n=bbs/Statistics/M.1617770863.A.74C.html
1F:推 andrew43: 推向量化,但對數學的基本功要求高不少。04/07 15:06
但是統計跟寫code本來就是需要大量數學的XDDDD 不要跟我說數學不好 可以把這兩個都學好(茶 所以有些現實還是要早一點讓學習者認清 數學不好,程式跟統計都很難到達top的水準 雖說不是每個人都需要到top的水準 但能有機會進步,誰不願意呢XDD
2F:推 joshddd: 推 不過有時候就是懶吧04/08 07:33
懶惰才更該把程式寫好(茶 技術債才不會越積越多,後面才不會改一個東西動輒得咎(攤手 再來就是這是寫code時間的花費與計算時間之間的取捨 你花多一點時間寫好的架構跟速度快的程式 跟花少少時間寫能動速度不快的程式 雖說都能得到一樣的結果 但是這套程式 你日後是不是還要再使用、會不會要修改它會有很大的不同 如果今日懶惰,之後又要新增很多功能的時候 後者保證讓你痛不欲生(笑 ※ 編輯: celestialgod (101.10.76.194 臺灣), 04/08/2021 09:35:20
3F:推 joshddd: 謝謝大大的建議 我會改進的 04/08 17:20
4F:→ DIDIMIN: 懶惰才是寫程式的原動力,只想寫出可以一鍵完成的程式 04/09 02:16







like.gif 您可能會有興趣的文章
icon.png[問題/行為] 貓晚上進房間會不會有憋尿問題
icon.pngRe: [閒聊] 選了錯誤的女孩成為魔法少女 XDDDDDDDDDD
icon.png[正妹] 瑞典 一張
icon.png[心得] EMS高領長版毛衣.墨小樓MC1002
icon.png[分享] 丹龍隔熱紙GE55+33+22
icon.png[問題] 清洗洗衣機
icon.png[尋物] 窗台下的空間
icon.png[閒聊] 双極の女神1 木魔爵
icon.png[售車] 新竹 1997 march 1297cc 白色 四門
icon.png[討論] 能從照片感受到攝影者心情嗎
icon.png[狂賀] 賀賀賀賀 賀!島村卯月!總選舉NO.1
icon.png[難過] 羨慕白皮膚的女生
icon.png閱讀文章
icon.png[黑特]
icon.png[問題] SBK S1安裝於安全帽位置
icon.png[分享] 舊woo100絕版開箱!!
icon.pngRe: [無言] 關於小包衛生紙
icon.png[開箱] E5-2683V3 RX480Strix 快睿C1 簡單測試
icon.png[心得] 蒼の海賊龍 地獄 執行者16PT
icon.png[售車] 1999年Virage iO 1.8EXi
icon.png[心得] 挑戰33 LV10 獅子座pt solo
icon.png[閒聊] 手把手教你不被桶之新手主購教學
icon.png[分享] Civic Type R 量產版官方照無預警流出
icon.png[售車] Golf 4 2.0 銀色 自排
icon.png[出售] Graco提籃汽座(有底座)2000元誠可議
icon.png[問題] 請問補牙材質掉了還能再補嗎?(台中半年內
icon.png[問題] 44th 單曲 生寫竟然都給重複的啊啊!
icon.png[心得] 華南紅卡/icash 核卡
icon.png[問題] 拔牙矯正這樣正常嗎
icon.png[贈送] 老莫高業 初業 102年版
icon.png[情報] 三大行動支付 本季掀戰火
icon.png[寶寶] 博客來Amos水蠟筆5/1特價五折
icon.pngRe: [心得] 新鮮人一些面試分享
icon.png[心得] 蒼の海賊龍 地獄 麒麟25PT
icon.pngRe: [閒聊] (君の名は。雷慎入) 君名二創漫畫翻譯
icon.pngRe: [閒聊] OGN中場影片:失蹤人口局 (英文字幕)
icon.png[問題] 台灣大哥大4G訊號差
icon.png[出售] [全國]全新千尋侘草LED燈, 水草

請輸入看板名稱,例如:e-shopping站內搜尋

TOP