読者です 読者をやめる 読者になる 読者になる

餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

Rは行列演算≿apply関数≿ループ演算

Rのプロビット回帰を用いたマイクロ・ベンチマークを書いたので、ソースコードを公開しておきます。

1. セットアップ

推定結果が常に同じになるように、セットアップします。

### セットアップ ###
set.seed(20120209)
obs <- 1000
crt <- runif(obs)
x1 <- runif(obs, 0, 1)
x2 <- runif(obs, 0, 2)
y <- ifelse(runif(obs) < pnorm(-1 + 0.7*x1 + 0.5*x2), 1, 0)
frml <- y ~ x1 + x2
# 説明変数を行列形式に変換
X <- model.matrix(frml)

2. 尤度関数の比較

対数尤度の計算関数において、行列演算、apply関数、ループ演算の3パターンを試します。

### 行列演算 ###
gc();gc()
system.time({
  r_nlm <- nlm(function(b){
    prd <- pnorm(X %*% matrix(b, ncol=1))
    # nlm()関数用に符号を逆にしておく
    -1*(sum(y*log(prd)) + sum*1
    })
    -1*(sum(y*log(prd)) + sum*2
    for(m in 1:nrow(X)){
      x <- X[m, ]
      p <- 0
      for(n in 1:length(b)){
        p <- p + b[n] * x[n]
      }
    prd[m] <- pnorm(p)
  }
  # nlm()関数用に符号を逆にしておく
  -1*(sum(y*log(prd)) + sum((1-y)*log(1-prd)))
}, c(0, 0, 0), hessian=TRUE)
})

3. MCMC(M-Hアルゴリズム)での実行

尤度関数の書き方でMCMCの実行時間は大きく異なる事になりますが、試した中で最速なものになります。

library(MCMCpack)
gc();gc()
system.time({
  post.samp <- MCMCmetrop1R(function(b, y, X){
    prd <- pnorm(X %*% matrix(b, ncol=1))
    sum(y*log(prd)) + sum( (1-y)*log(1-prd) )
  }, theta.init=c(0,0,0), X=X, y=y, logfun=TRUE)
})

*1:1-y)*log(1-prd))) }, c(0, 0, 0), hessian=TRUE) }) ### apply() ### gc();gc() system.time({   r_nlm <- nlm(function(b){     prd <- apply(X, 1, function(n){       r <- pnorm(sum(n*b

*2:1-y)*log(1-prd))) }, c(0, 0, 0), hessian=TRUE) }) ### ループ演算 ### gc();gc() system.time({   r_nlm <- nlm(function(b){     prd <- numeric(nrow(X