作者celestialgod (天)
看板R_Language
标题Re: [问题] 单一y值针对每两个栏位去计算lm
时间Fri Oct 7 21:31:49 2016
※ 引述《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://webptt.com/cn.aspx?n=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
1F：推 cywhale: 收起来学，谢谢分享^^~~ 10/07 21:40
2F：推 andrew43: XD 10/07 21:42
3F：推 Edster: 这真的要学起来。 10/07 21:49
多提供一个好了，LOOCV in RcppParallel
http://pastebin.com/4JK6VSd7
※ 编辑: celestialgod (36.233.137.163), 10/07/2016 21:51:52
4F：推 yu8861213: 版版好帅^Q^Q^Q^Q^ 10/07 21:59