重回帰はなぜ分散処理できる? - 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()より処理速度が遅いです。(笑)
あと、一台のマシンですべての処理をしているので、
テーブルを細かく分ければ分けるほど、処理時間は長くなります(笑)