雙迴路中的線性模型？

當n增大時，R中double循環的速度非常慢。有沒有辦法通過for循環提高速度？雙迴路中的線性模型？

set.seed(1) 
    n=1000 

    y=rnorm(n) 
    x1=rnorm(n) 
    x2=rnorm(n) 

    lm.ft=function(y,x1,x2) 
     lm.fit(cbind(1,x1.bar,x2.bar), y)$coef 

    res=array(,dim=c(1,3,n,n)) 
    for(i in 1:n) 
     for(j in 1:n){ 
     x1.bar=x1-x1[i] 
     x2.bar=x2-x2[j] 
     res[,,i,j]=lm.ft(y,x1.bar,x2.bar) 
     }

來源

2013-12-19 user1690124

你不可能期望它在'n'中的增長幅度小於二次方。這個例子將需要3,000,000個單獨的調用到'lm.fit'。你應該解釋你正在做什麼，（或被要求做），並請指出這是你的家庭作業問題。 –

@IShouldBuyABoat我只給lm.fit打了100萬個電話，但這肯定夠用了。 OP沒有意識到98％以上的時間都花在了'lm.fit'裏面。 –

同意。我正在尋找第二個指數來暗淡（）和相乘。如果我能想出實際的目標，我想這可能是通過對單個lm.fit進行簡單的增量調整來解決的。 –

只給你一個完整的答案：除了在你的代碼中的一些古怪（如在使用x1.bar和x2.bar內lm.ft而不是x1和x2），我分析說：你到底想達到什麼????

如果我對自己的代碼運行此：

Rprof("profile1.out") 
for(i in 1:n) 
    for(j in 1:n){ 
    x1.bar=x1-x1[i] 
    x2.bar=x2-x2[j] 
    res[,,i,j]=lm.ft(y,x1.bar,x2.bar) 
    } 
Rprof(NULL) 
summaryRprof("profile1.out")

我得到以下有趣的畫面：你只是擬合模型時

> summaryRprof("profile1.out") 
$by.self 
       self.time self.pct total.time total.pct 
".Call"    0.96 22.86  0.96  22.86 
"lm.fit"    0.92 21.90  4.08  97.14 
... 
"cbind"    0.22  5.24  0.22  5.24 
... 

$by.total 
       total.time total.pct self.time self.pct 
"lm.ft"    4.12  98.10  0.04  0.95 
"lm.fit"    4.08  97.14  0.92 21.90 
... 
"cbind"    0.22  5.24  0.22  5.24 
...

98％。循環速度並不慢，事實上，您正在嘗試安裝100萬個模型，這讓您等待。你真的不得不重新思考你的問題。

如果這真的是你想要做的，那麼優化你的功能將涉及到擺脫lm.fit的開銷和矢量化減法。節省約50％。

lm.ft=function(y,x1,x2) 
    .Call(stats:::C_Cdqrls, cbind(1,x1,x2), y, tol=1e-7)$coef 

x1.bar <- outer(x1,x1,`-`) 
x2.bar <- outer(x2,x2,`-`) 
for(i in 1:n) 
    for(j in 1:n){ 
    res[,,i,j]=lm.ft(y,x1.bar[,i],x2.bar[,j]) 
    }

來源

2013-12-19 17:43:24

生存套餐中是否存在coxph的優化方式？代碼是coxph（Surv（t，delta）〜z1 + z2）$係數 – user1690124

這回到了你想要做什麼的問題。如果你提供了更多的上下文，人們可能會給出解決方案，以避免需要適應100萬個模型... –

我的意思是像.Call的代碼（stats ::: C_Cdqrls，cbind（1，x1，x2），y ，tol = 1e-7）$ coef。 – user1690124

如果你想要做一些瘋狂的樣子，你應該使用RCPP：

library(RcppEigen) 
library(inline) 

incl <- ' 
using Eigen::LLT; 
using Eigen::Lower; 
using Eigen::Map; 
using Eigen::MatrixXd; 
using Eigen::MatrixXi; 
using Eigen::Upper; 
using Eigen::VectorXd; 
using Eigen::Vector3d; 
typedef Map<MatrixXd> MapMatd; 
typedef Map<MatrixXi> MapMati; 
typedef Map<VectorXd> MapVecd; 
inline MatrixXd AtA(const MatrixXd& A) { 
    int n(A.cols()); 
    return MatrixXd(n,n).setZero().selfadjointView<Lower>().rankUpdate(A.adjoint()); 
} 
' 

body <- ' 
const MapMatd  X(as<MapMatd>(XX)); 
const MapVecd  y(as<MapVecd>(yy)); 
const int   n(X.rows()), m(X.cols()); 
LLT<MatrixXd>  llt; 
MatrixXd    Res(n*n,m), Xbar(n,m); 
Vector3d    betahat; 
for (int i = 0; i < n; ++i) { 
for (int j = 0; j < n; ++j) { 
    Xbar=X; 
    for (int k = 0; k < n; ++k) { 
    Xbar(k,1) -= X(i,1); 
    Xbar(k,2) -= X(j,2); 
    }; 
    llt=AtA(Xbar); 
    betahat =llt.solve(Xbar.adjoint() * y); 
    Res.row(i*n+j) = betahat; 
}; 
}; 
return    wrap(Res); 
' 

crazyLm <- cxxfunction(signature(XX = "matrix", yy = "numeric"), 
          body, "RcppEigen", incl) 

set.seed(1) 
n=4 

y=rnorm(n) 
x1=rnorm(n) 
x2=rnorm(n) 

lm.ft=function(y,x1,x2) lm.fit(cbind(1,x1.bar,x2.bar), y)$coef 

res=array(,dim=c(3,n,n)) 
for(i in 1:n) 
    for(j in 1:n){ 
    x1.bar=x1-x1[i] 
    x2.bar=x2-x2[j] 
    res[,i,j]=lm.ft(y,x1.bar,x2.bar) 
    } 

res2 <- aperm(array(t(crazyLm(cbind(1,x1,x2), y)), dim=c(3,n,n)), c(1,3,2)) 
all.equal(res, res2) 
#[1] TRUE 

system.time({ 
set.seed(1) 
n=1000 

y=rnorm(n) 
x1=rnorm(n) 
x2=rnorm(n) 
res <- aperm(array(t(crazyLm(cbind(1,x1,x2), y)), dim=c(3,n,n)), c(1,3,2)) 
}) 

# User  System  elapsed 
#36.130  0.033  36.158

這可以讓你適應百萬車型在不到一分鐘。但是，我沒有看到一個用例。

來源

2013-12-19 21:20:27 Roland

雙迴路中的線性模型？

回答

相關問題