作者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/m.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