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
