重回帰はなぜ分散処理できる？ - qz70224の統計学メモ
のロジスティック回帰版です。
根本の考え方は基本的に重回帰と同じなので、詳しくはそちらをご覧ください。
段取りは以下の通りです。
ニュートン法については、以下の資料を参考にさせていただきました。
http://www.chokkan.org/publication/survey/prml_chapter4_discriminative_slides.pdf
以下、プログラムです。
サンプルデータ作成
library(MASS)
mu<-c(0,0,0,0)
Sigma<-matrix(NA,4,4)
diag(Sigma)<-1
Sigma[lower.tri(Sigma)]<-c(0.4,0.3,0.1,0.2,0.3,0.6)
Sigma[upper.tri(Sigma)]<-c(0.4,0.3,0.1,0.2,0.3,0.6)
data<-mvrnorm(n=1000000,mu=mu,Sigma=Sigma,empirical=F)
for(m in 1:nrow(data)){
if(pnorm(data[m,4])>0.5){
data[m,4]=1
}else{
data[m,4]=0
}
}
lst<-split(as.data.frame(data),1:10)計算
#lst=リスト形式のデータ,num=従属変数になる列番号
mpLR<-function(lst,num){
#デザイン行列を作成する（対象列番号を除外して切片をつける）
dzmatlst<-split(rep(NA,length(lst)),1:length(lst),1)
for(k in 1:length(dzmatlst)){
dzmatlst[[k]] <- as.matrix( cbind(1,lst[[k]][,-num]) )
}
#リスト毎に重回帰を実施して初期値決定
#リスト数×パラメータ数の行列を作成
wlm<-matrix(NA,ncol(dzmatlst[[1]]),length(lst))
#リスト毎に重回帰実施し重みを行列に格納
for(a in 1:length(lst)){
wlm[,a]<-solve( t(dzmatlst[[a]])%*%dzmatlst[[a]] ) %*% t(dzmatlst[[a]]) %*% lst[[a]][,num]
}
#他に特にいい方法が思いつかなかったので、すべての重回帰のパラメータ推定値を平均して初期値にした
w<-apply(wlm,1,mean)
#以下よりループ処理
repeat{
yhat <- split(rep(NA,length(lst)),1:length(lst),1)
dfb <- matrix(NA,ncol(dzmatlst[[1]]),length(lst))
hess <- array(NA,dim=c(ncol(dzmatlst[[1]]),ncol(dzmatlst[[1]]),length(lst)))
for(a in 1:length(lst)){
#リスト毎に予測値定義
yhat[[a]] <- 1 / ( 1 + exp(-1 * (dzmatlst[[a]] %*% w) ) )
#リスト毎に勾配
dfb[,a] <- t(dzmatlst[[a]]) %*% (yhat[[a]]-lst[[a]][,num])
#リスト毎にヘッセ行列
R<-yhat[[a]]*(1-yhat[[a]])
tXR<-matrix(NA,ncol(dzmatlst[[a]]),nrow(dzmatlst[[a]]))
for(k in 1:nrow(dzmatlst[[a]])){
tXR[,k]<-dzmatlst[[a]][k,]*R[k]
}
hess[,,a]<-tXR%*%dzmatlst[[a]]
}
#勾配、ヘッセ行列の和を求める
dfb<-rowSums(dfb)
hess<-apply(hess,c(1,2),sum)
#重みの更新
w1 <- w-solve(hess)%*%dfb
#収束判定
if(max(w1-w)>0.0001){
w<-w1
}else{
break
}
}
se_w<-sqrt(diag(solve(hess)))
tval<-w/se_w
for(l in 1:length(lst)){
if(l==1){
loglik<-sum( log(yhat[[l]])*lst[[l]][,num] + log(1-yhat[[l]])*(1-lst[[l]][,num]) )
}else{
loglik<-loglik + sum( log(yhat[[l]])*lst[[l]][,num] + log(1-yhat[[l]])*(1-lst[[l]][,num]) )
}
}
AIC<--2*(loglik-ncol(dzmatlst[[1]]))
#出力結果
out<-list(t(w),se_w,t(tval),-2*loglik,AIC)
names(out)<-c("w","se_w","tval","-2loglik","AIC")
return(out)
}試してみる
start<-proc.time()
res<-mpLR(lst,4)
mpLRtime<-proc.time()-start
> res
$w
[,1] [,2] [,3] [,4]
[1,] 0.0004474997 -0.4053877 0.5541994 1.335781
$se_w
[1] 0.002353378 0.002715095 0.002707285 0.003178775
$tval
[,1] [,2] [,3] [,4]
[1,] 0.1901521 -149.3088 204.7067 420.2187
$`-2loglik`
[1] 1073944
$AIC
[1] 1073952
> mpLRtime
ユーザ システム 経過
34.69 0.06 35.12
start<-proc.time()
LR<-glm(data[,4]~data[,1]+data[,2]+data[,3],
family=binomial(link="logit"))
summary(LR)
LRtime<-proc.time()-start
> summary(LR)
Call:
glm(formula = data[, 4] ~ data[, 1] + data[, 2] + data[, 3],
family = binomial(link = "logit"))
Deviance Residuals:
Min 1Q Median 3Q Max
-3.1390 -0.8898 0.1131 0.8901 3.0114
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.0004475 0.0023534 0.19 0.849
data[, 1] -0.4053878 0.0027151 -149.31 <2e-16 ***
data[, 2] 0.5541995 0.0027072 204.71 <2e-16 ***
data[, 3] 1.3357808 0.0031787 420.23 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1386294 on 999999 degrees of freedom
Residual deviance: 1073944 on 999996 degrees of freedom
AIC: 1073952
Number of Fisher Scoring iterations: 4パラメータ推定値は一致しました。
glm()より処理速度が遅いです。（笑）
あと、一台のマシンですべての処理をしているので、
テーブルを細かく分ければ分けるほど、処理時間は長くなります（笑）
