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/cn.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灯, 水草

请输入看板名称,例如:Gossiping站内搜寻

TOP