Re: [問題] 多筆資料運算Pearson相關系數

看板R_Language作者celestialgod (攸藍)時間9年前 (2015/06/11 12:49)推噓2(2推 0噓 3→)

留言5則, 4人參與討論串3/3 (看更多)

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://www.ptt.cc/bbs/R_Language/M.1433998170.A.218.html