Re: [問題] 單一y值針對每兩個欄位去計算lm

看板R_Language作者celestialgod (天)時間9年前 (2016/10/07 21:31)推噓4(4推 0噓 0→)

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

※ 引述《girl5566 (5566520)》之銘言： : 文章分類提示: : - 問題: 當你想要問問題時，請使用這個類別 : [問題類型]: : 程式諮詢(我想用R 做某件事情，但是我不知道要怎麼用R 寫出來) : [軟體熟悉度]: : 新手(沒寫過程式，R 是我的第一次) : [問題敘述]: : 想嘗試將所有的X組合跟Y進行lm計算，目前不太清楚何種方式較有效率 : 是否應該要用平行處理!? : [程式範例]: : fun <- function(x) coef(lm(paste0("y ~ ", paste(x, collapse="+")), : data=iris))[2] : combn(names(iris[1:4]), 2, fun) : [環境敘述]: : Win7 64 bit , R 3.2.2 : [關鍵字]: 最快的方法： library(Rcpp) library(RcppArmadillo) library(RcppParallel) sourceCpp(code = " // [[Rcpp::depends(RcppArmadillo)]] #include <RcppArmadillo.h> using namespace Rcpp; using namespace arma; // [[Rcpp::depends(RcppParallel)]] #include <RcppParallel.h> using namespace RcppParallel; struct coef_compute: public Worker { const mat X; const vec& y; const umat& combnMat; mat& coefMat; coef_compute(const mat& X, const vec& y, const umat& combnMat,mat& coefMat): X(X), y(y), combnMat(combnMat), coefMat(coefMat) {} void operator()(std::size_t begin, std::size_t end) { for (uword i = begin; i < end; i++) coefMat.row(i) = solve(join_rows(ones<vec>(y.n_elem), X.cols(combnMat.col(i))), y).t(); } }; // [[Rcpp::export]] List fastCoef_allCombns(arma::vec y, arma::mat x, int k) { Rcpp::Function combn_rcpp(\"combn\"); umat combnMat = conv_to<umat>::from(as<mat>(combn_rcpp(x.n_cols, k))) - 1; mat coefMat = zeros<mat>(combnMat.n_cols, k + 1); coef_compute coefResults(x, y, combnMat, coefMat); parallelFor(0, combnMat.n_cols, coefResults); return List::create(combnMat.t() + 1, coefMat); }") n <- 30 X <- matrix(rnorm(120), n) coefs <- c(2, 5, -7, 8, 0.5) y <- as.vector(cbind(1, X) %*% coefs) + rnorm(n, 0, 5) outRes <- lapply(1:ncol(X), function(i){ fastCoef_allCombns(y, X, i) }) outRes[[2]] # # 哪些X的組合 # [[1]] # [,1] [,2] # [1,] 1 2 # X1 跟 X2的組合 # [2,] 1 3 # [3,] 1 4 # [4,] 2 3 # [5,] 2 4 # [6,] 3 4 # # # 每一個組合的X，其係數，BY列去看 # [[2]] # [,1] [,2] [,3] # [1,] 0.4210033 6.445250 -7.5248928 # [2,] 0.7895878 7.378273 8.7397441 # [3,] 0.6180885 6.864235 0.9617974 # [4,] 2.6213473 -6.677914 6.3579295 # [5,] 2.3471788 -7.834296 1.4096750 # [6,] 2.9802756 8.212602 0.9764216 # benchmark n <- 3000; p <- 10 X <- matrix(rnorm(n*p), n) coefs <- sample(-10:10, p+1, TRUE) y <- as.vector(cbind(1, X) %*% coefs) + rnorm(n, 0, 5) st <- proc.time() outRes <- lapply(1:ncol(X), function(i){ fastCoef_allCombns(y, X, i) }) proc.time() - st # user system elapsed # 0.14 0.08 0.07 st <- proc.time() outRes2 <- lapply(1:ncol(X), function(i){ combnMat <- combn(1:ncol(X), i) return(list(t(combnMat), t(apply(combnMat, 2, function(v){ coef(lm(y ~ X[,v])) })))) }) proc.time() - st # user system elapsed # 3.57 0.19 3.56 快沒很多啦，快50倍而已 -- R資料整理套件系列文： magrittr #1LhSWhpH (R_Language) http://tinyurl.com/j3ql84c data.table #1LhW7Tvj (R_Language) http://tinyurl.com/hr77hrn dplyr(上) #1LhpJCfB (R_Language) http://tinyurl.com/jtg4hau dplyr(下) #1Lhw8b-s (R_Language) tidyr #1Liqls1R (R_Language) http://tinyurl.com/jq3o2g3 -- ※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 36.233.137.163 ※ 文章網址: https://www.ptt.cc/bbs/R_Language/M.1475847112.A.D6B.html ※ 編輯: celestialgod (36.233.137.163), 10/07/2016 21:38:23 ※ 編輯: celestialgod (36.233.137.163), 10/07/2016 21:39:57