4

"lm" を使用せずに R で通常の最小二乗 ( OLS ) 推定値を計算したいのですが、これにはいくつかの理由があります。まず、「lm」は、私の場合はデータ サイズが問題であることを考慮して、必要のない多くのもの (適合値など) も計算します。第二に、OLS を別の言語で実装する前に R で自分自身を実装できるようにしたい (たとえば C で GSL を使用する)。

ご存知かもしれませんが、モデルは次のとおりです。Y=Xb+E; E〜N（0、シグマ^ 2）で。以下で詳しく説明するように、b は平均 (b0) と別の係数 (b1) の 2 つのパラメーターを持つベクトルです。最後に、実行する線形回帰ごとに、b1 (効果サイズ) の推定値、その標準誤差、sigma^2 (残差分散)、および R^2 (決定係数) の推定値が必要です。

これが私のデータです。N 個のサンプルがあります (例: 個人、N~=100)。サンプルごとに、Y 出力 (応答変数、Y~=10^3)​​ と X ポイント (説明変数、X~=10^6) があります。Y出力を別々に扱いたい、つまり. Y 線形回帰を開始したい: 出力 1 用に 1 つ、出力 2 用に 1 つなど。それから...最終的にポイントXに。

これは、 「lm」の速度と、行列代数を介した OLS 推定の計算速度をチェックするための私のR コードです。

まず、ダミー データをシミュレートします。

nb.samples <-  10  # N
nb.points <- 1000  # X
x <- matrix(data=replicate(nb.samples,sample(x=0:2,size=nb.points, replace=T)),
            nrow=nb.points, ncol=nb.samples, byrow=F,
            dimnames=list(points=paste("p",seq(1,nb.points),sep=""),
              samples=paste("s",seq(1,nb.samples),sep="")))
nb.outputs <- 10  # Y
y <- matrix(data=replicate(nb.outputs,rnorm(nb.samples)),
            nrow=nb.samples, ncol=nb.outputs, byrow=T,
            dimnames=list(samples=paste("s",seq(1,nb.samples),sep=""),
              outputs=paste("out",seq(1,nb.outputs),sep="")))

すぐ下で使用される私自身の関数は次のとおりです。

GetResFromCustomLinReg <- function(Y, xi){ # both Y and xi are N-dim vectors
  n <- length(Y)
  X <- cbind(rep(1,n), xi)  #
  p <- 1      # nb of explanatory variables, besides the mean
  r <- p + 1  # rank of X: nb of indepdt explanatory variables
  inv.XtX <- solve(t(X) %*% X)
  beta.hat <- inv.XtX %*% t(X) %*% Y
  Y.hat <- X %*% beta.hat
  E.hat <- Y - Y.hat
  E2.hat <- (t(E.hat) %*% E.hat)
  sigma2.hat <- (E2.hat / (n - r))[1,1]
  var.covar.beta.hat <- sigma2.hat * inv.XtX
  se.beta.hat <- t(t(sqrt(diag(var.covar.beta.hat))))
  Y.bar <- mean(Y)
  R2 <- 1 - (E2.hat) / (t(Y-Y.bar) %*% (Y-Y.bar))
  return(c(beta.hat[2], se.beta.hat[2], sigma2.hat, R2))
}

組み込みの「lm」を使用したコードは次のとおりです。

res.bi.all <- apply(x, 1, function(xi){lm(y ~ xi)})

これが私のカスタム OLS コードです。

res.cm.all <- apply(x, 1, function(xi){apply(y, 2, GetResFromCustomLinReg, xi)})

上記の値でこの例を実行すると、次のようになります。

> system.time( res.bi.all <- apply(x, 1, function(xi){lm(y ~ xi)}) )
   user  system elapsed
  2.526   0.000   2.528
> system.time( res.cm.all <- apply(x, 1, function(xi){apply(y, 2, GetResFromCustomLinReg, xi)}) )
   user  system elapsed
  4.561   0.000   4.561

(そして、当然、N、X、Y を増やすと悪化します。)

もちろん、"lm" には応答行列 (y~xi) の各列を個別に "自動的に" フィッティングするという優れた特性がありますが、Y 回 (yi~xi ごとに) "適用" する必要があります。しかし、これが私のコードが遅い唯一の理由ですか? これを改善する方法を知っている人はいますか？

(非常に長い質問で申し訳ありませんが、最小限でありながら包括的な例を提供しようとしました。)

> sessionInfo()
R version 2.12.2 (2011-02-25)
Platform: x86_64-redhat-linux-gnu (64-bit)
4

2 に答える 2

7

CRAN のRcppArmadilloパッケージのfastLm()関数を見てください。これに先行するRcppGSLにも同様のものがありますが、おそらくArmadilloベースのソリューションが必要です。古いプレゼンテーション (R を使用した HPC) には、速度の向上を示すスライドがいくつかあります。fastLm()

また、縮退したモデル行列で問題になる可能性がある X'X の直接の逆よりも優れた「ピボット」アプローチに関するヘルプ ページのヒントにも注意してください。

于 2011-06-04T04:24:06.233 に答える
4

Marek のコメントに従って、組み込み関数「lm」と「lm.fit」、私自身の関数、パッケージ RcppArmadillo の「fastLm」と「fastLmPure」を比較した結果を以下に示します。

> system.time( res1 <- apply(x, 1, function(xi){lm(y ~ xi)}) )
   user  system elapsed
  2.859   0.005   2.865
> system.time( res2 <- apply(x, 1, function(xi){apply(y, 2, GetResFromCustomLinReg, xi)}) )
   user  system elapsed
  4.620   0.004   4.626
> system.time( res3 <- apply(x, 1, function(xi){lm.fit(x=cbind(1,xi), y=y)}) )
   user  system elapsed
  0.454   0.004   0.458
> system.time( res4 <- apply(x, 1, function(xi){apply(y, 2, fastLm, x=cbind(1,xi))}) )
   user  system elapsed
  2.279   0.005   2.283
> system.time( res5 <- apply(x, 1, function(xi){apply(y, 2, fastLmPure, cbind(1,xi))}) )
   user  system elapsed
  1.053   0.003   1.056

ただし、これらの数値を比較する場合は注意が必要です。違いは、実装の違いだけでなく、結果が効果的に計算されることにも起因します。

> names(res1$p1)
 [1] "coefficients"  "residuals"     "effects"       "rank"        
 [5] "fitted.values" "assign"        "qr"            "df.residual" 
 [9] "xlevels"       "call"          "terms"         "model"       
> # res2 (from my own custom function) returns the estimate of beta, its standard error, the estimate of sigma and the R^2
> names(res3$p1)
[1] "coefficients"  "residuals"     "effects"       "rank"        
[5] "fitted.values" "assign"        "qr"            "df.residual" 
> names(res4$p1$out1)
[1] "coefficients"  "stderr"        "df"            "fitted.values"
[5] "residuals"     "call"        
> names(res5$p1$out1)
[1] "coefficients" "stderr"       "df"

たとえば、「lm」よりも「lm.fit」を使用することを好むかもしれませんが、R^2 が必要な場合は、自分で計算する必要があります。同上、「lm」よりも「fastLm」を使用したい場合がありますが、シグマの推定値が必要な場合は、自分で計算する必要があります。また、カスタム R 関数を使用してそのような計算を行うのは、あまり効率的ではない可能性があります (「lm」で行われることと比較してください)。

これらすべてに照らして、当面は「lm」を使用し続けますが、実際、「fastLm」に関するDirkのコメントは良いアドバイスです（他の人にとって興味深いはずなので、彼の答えを選んだのはそのためです）。

于 2011-06-27T16:24:37.387 に答える