作者celestialgod (攸蓝)
看板R_Language
标题Re: [问题] 多笔资料运算Pearson相关系数
时间Thu Jun 11 12:49:26 2015
2018/06/14更新，全面修改程式，增加一些范例，url：
https://pastebin.com/ssqrKBcG
library(magrittr)
library(data.table)
library(reshape2)
library(tidyr)
library(dplyr)
## data generation
N = 5e5
dat = data.table(paste("sample", 1:N), matrix(rnorm(N*60),, 60))
dat %<>% setnames(c("Name", paste0(rep(c("X", "Y"),,, 30), rep(1:30,2)))) %>%
tbl_dt(FALSE)
## for loop
st = proc.time()
corCoef = vector('numeric', N)
for (i in 1:N){
corCoef[i] = cor(as.numeric(dat[i, 2:31, with=FALSE]),
as.numeric(dat[i, 32:61, with=FALSE]))
}
output = data.table(Name = dat$Name, cor = corCoef)
proc.time() - st
# user system elapsed
# 1345.32 0.80 1351.27
## melt
st = proc.time()
dat2 = melt(dat, id = 1) %>% mutate(G = as.integer(grepl("X", variable)))
output2 = dat2 %>% group_by(Name) %>%
summarise(cor = cor(value[G==1], value[G==0]))
proc.time() - st
# user system elapsed
# 30.81 0.33 31.26
all.equal(output$cor, output2$cor)
# [1] TRUE
用对方法会快上很多！
两个方法大概差了43倍，学会把资料做适当的变换，
并加以分组，用分组的方式去计算会比较快速而有效率！
这部分也包含你资料一开始整理的方式不是很恰当。
版友提供：
library(matrixStats)
dat3 = as.matrix(dat[,2:61, with=FALSE])
st = proc.time()
m_X <- dat3[,1:30]
m_Y <- dat3[,31:60]
mu_X <- rowMeans(m_X)
mu_Y <- rowMeans(m_Y)
sigma_X <- rowSds(m_X)
sigma_Y <- rowSds(m_Y)
myCorr <- rowSums((m_X-mu_X)*(m_Y-mu_Y))/((sigma_X*sigma_Y)*(30-1))
proc.time() - st
# user system elapsed
# 0.50 0.06 0.56
all.equal(output2$cor, myCorr)
# [1] TRUE
## Rcpp
library(Rcpp)
library(RcppArmadillo)
sourceCpp(code = '
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
NumericVector fastCor(NumericMatrix Xr, NumericMatrix Yr) {
uword k = Xr.ncol();
mat X(Xr.begin(), Xr.nrow(), k, false);
mat Y(Yr.begin(), Xr.nrow(), k, false);
colvec output(Xr.nrow());
X.each_col() -= mean(X, 1);
X.each_col() /= stddev(X, 0, 1);
Y.each_col() -= mean(Y, 1);
Y.each_col() /= stddev(Y, 0, 1);
output = sum(X % Y, 1) / ((double) k - 1);
return wrap(output);
}')
st = proc.time()
output3 = fastCor(dat3[,1:30], dat3[,31:60])
proc.time() - st
# user system elapsed
# 0.44 0.01 0.45
all.equal(output2$cor, as.vector(output3))
# [1] TRUE
※ 引述《Shadowy (+ing》之铭言：
: 各位先进好：
: 我有50万笔资料，每一笔资料有30组(X,Y)的数据，
: 想要针对每一笔的X,Y运算Pearson相关系数，
: 资料格式，如下：
: Name X1 X2 ... X30 Y1 Y2 ... Y30
: .
: .
: .
: .
: 共50万笔
: 欲输出格式为：
: (Name) (Pearson's cor)
: 因为没有太多的程式撰写经验，
: 目前的想法是：
: 先抓取每一列1~30个值为X向量，31~60个值为Y向量，
: 进行cor(X, Y, use="complete", method="pearson")运算Pearson相关系数，
: 再利用回圈运算50万笔资料。
: 请问先进，我应该如何开始撰写这样的语法呢？
: 还是我应该改变汇入资料的格式呢？
: 再麻烦各位先进指教！
: 谢谢大家~
--
※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 36.225.240.107
※ 文章网址: https://webptt.com/cn.aspx?n=bbs/R_Language/M.1433998170.A.218.html
1F：→ Shadowy: 谢谢c大解惑与建议，我会再努力，希望有朝一日我也能像各 06/11 16:12
2F：→ Shadowy: 位一样。 06/11 16:12
3F：推 cywhale: great to evaluate by original formula >> cor() ^^ 06/11 22:30
4F：→ andrew43: 0.5秒…太神了 06/12 19:12
※ 编辑: celestialgod (36.225.212.153), 06/12/2015 20:45:07
5F：推 Edster: 蛮值得参考的，谢谢。 06/14 10:09
※ 编辑: celestialgod (111.246.26.103), 06/14/2018 00:51:00