データの読み込み
#以下のパッケージをインストールしておいてください．
library(ISLR)
library(dplyr)
パッケージISLR内のDefaultデータを用いる．説明変数をカード借入金額(balance)，結果変数を債務不履行(default)とし，結果変数はダミー変数(0 or 1)にしておく．データ処理を行うために，パッケージdplyrも読み込んでおく．
defa <- Default[,c("default","balance")] #ISLR内のDefaultデータを使う
defa <- mutate(defa,default=as.numeric(default=="Yes")) #defalt変数をダミー化する
head(defa)
| default | balance | |
|---|---|---|
| <dbl> | <dbl> | |
| 1 | 0 | 729.5265 |
| 2 | 0 | 817.1804 |
| 3 | 0 | 1073.5492 |
| 4 | 0 | 529.2506 |
| 5 | 0 | 785.6559 |
| 6 | 0 | 919.5885 |
散布図を描いてみる
plot(x=NULL,y=NULL,xlim=c(0,max(defa[,2]+1)),ylim=c(0,1),xlab="Balance",ylab="Default")
points(defa[which(defa$default==1),2],defa[which(defa$default==1),1],col="blue")
points(defa[which(defa$default==0),2],defa[which(defa$default==0),1],col="red")
トレーニングデータとテストデータ
defa.tr <- defa[1:9000,]
defa.te <- defa[9001:10000,]
まずは線形回帰モデルを当てはめてみる
当てはめ後，クロス集計表を作成する．
reg.defa <- lm(default~balance,data=defa.tr) #トレーニングデータで学習
pre <- predict(reg.defa,defa.te) #テストデータに対して予測
table(defa.te[,1], as.numeric(pre > 0.5)) #クロス集計表の作成
0
0 964
1 36
plot(x=NULL,y=NULL,xlim=c(0,max(defa[,2]+1)),ylim=c(0,1),xlab="Balance",ylab="Default")
points(defa[which(defa$default==1),2],defa[which(defa$default==1),1],col="blue")
points(defa[which(defa$default==0),2],defa[which(defa$default==0),1],col="red")
abline(coef(reg.defa)[1],coef(reg.defa)[2],lwd=2)
ロジスティック回帰の実行
まずはRの関数に頼る方法．
log.defa <- glm(default~balance,data=defa.tr,family="binomial") #ロジスティック回帰モデルの当てはめ
summary(log.defa)
pre2 <- predict(log.defa,defa.te,type="response") #typeをresponseに指定することで確率を返す
table(defa.te[,1],as.numeric(pre2 > 0.5)) #クロス集計表の作成
Call:
glm(formula = default ~ balance, family = "binomial", data = defa.tr)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.2897 -0.1422 -0.0557 -0.0208 3.7175
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -10.817741 0.389584 -27.77 <2e-16 ***
balance 0.005596 0.000237 23.61 <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: 2610.4 on 8999 degrees of freedom
Residual deviance: 1402.1 on 8998 degrees of freedom
AIC: 1406.1
Number of Fisher Scoring iterations: 8
0 1
0 960 4
1 27 9
plot(x=NULL,y=NULL,xlim=c(0,max(defa[,2]+1)),ylim=c(0,1),xlab="Balance",ylab="Default")
points(defa.tr[which(defa.tr$default==1),2],defa.tr[which(defa.tr$default==1),1])
points(defa.tr[which(defa.tr$default==0),2],defa.tr[which(defa.tr$default==0),1])
#テストデータもプロット
points(defa.te[which(defa.te$default==1),2],defa.te[which(defa.te$default==1),1],col="blue",pch=3)
points(defa.te[which(defa.te$default==0),2],defa.te[which(defa.te$default==0),1],col="red",pch=3)
par(new=T)
#ロジスティック回帰による予測曲線
curve(1/(1+exp(-coef(log.defa)[1] - coef(log.defa)[2]*x)),xlim=c(0,max(defa[,2]+1)),ylim=c(0,1),xlab="",ylab="",lwd=2)
abline(a=0.5,b=0)
最尤推定によるパラメータの推定
ニュートン法で最適化を行う．
ff <- function(theta){1/(1+exp(-theta))}
x <- defa.tr[,2]; X <- cbind(1,x); y <- defa.tr[,1]
beta <- c(0,0); judge <- 1000
for(k in 1:100){ #高々何回繰り返すかを指定
if(judge >= 0.0001){ #値の更新がほとんど行われなくなったら終了
old <- beta
p <- ff(X%*%old)
H <- solve(t(X) %*% diag(c(p*(1-p))) %*% X)
beta <- old - H %*% t(X) %*% (p-y)
judge <- sum(abs(beta-old))
} else break
}
beta
| -10.817741116 | |
| x | 0.005595708 |
Exercise : 上記の方法だと$H$内のdiagを取る部分の行列サイズが大きく，計算に時間がかかってしまう．この点に注意しながら$H$を効率的に計算し，上記のコードを修正せよ．
Spotifyデータの解析をしながら説明する．
dtrain <- read.csv("disneylist_train.csv", header=T) #トレーニングデータの読み込み
title.train <- dtrain$track.name #曲名を保持しておく
dtrain <- dtrain[,-c(1,17,18,19)] #不要な変数を削除
dtrain <- mutate(dtrain,eval=as.numeric(eval>=5)) #evalが5以上を1, それ以外を0に変更
#テストデータに対しても同様の処理を行う
dtest <- read.csv("disneylist_test.csv", header=T)
title.test <- dtest$track.name
dtest <- dtest[,-c(1,17,18,19)]
dtest <- mutate(dtest,eval=as.numeric(eval>=5))
ロジスティック回帰の実行．glm関数のoptionにfamily=binomialを指定することで可能．もちろん，最尤推定で直接推定することもできる．
result <- glm(eval~danceability+energy+key,data=dtrain,family="binomial")
summary(result)
Call:
glm(formula = eval ~ danceability + energy + key, family = "binomial",
data = dtrain)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.1369 -0.7279 -0.6261 -0.4604 2.0858
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.21537 0.64714 -1.878 0.0604 .
danceability -2.15179 1.24558 -1.728 0.0841 .
energy 2.03711 0.95316 2.137 0.0326 *
key -0.07547 0.04794 -1.574 0.1154
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 241.83 on 231 degrees of freedom
Residual deviance: 234.60 on 228 degrees of freedom
AIC: 242.6
Number of Fisher Scoring iterations: 4
result <- glm(eval~.,data=dtrain,family="binomial") #全ての説明変数を使いたい場合
単変量の場合と同様，トレーニングデータで学習し，テストデータの中から曲を推薦してみる．
pre2 <- predict(result,dtest,type="response")
title.test[which(pre2 > 0.5)]
クロス集計表の作成．
table(dtest[,1],as.numeric(pre2 > 0.5))
0 1
0 22 2
1 12 0