ロジスティック回帰を導出から実装までみっちり途中式解説


#線形単回帰知ったボク

ふーん
最小二乗法便利じゃん


#ロジスティック回帰知ったボク

ふーん・・・
？？？なんでそれで係数求まるの？？？

ということで解説しながらやっていきます

#ロジスティック回帰についての疑問点5stepで。

##前提条件

ある調査によって男女の身長を求めた。
身長データの平均は160であった。
2行取り出した以下のデータを使いロジスティック回帰を行ってみよ。
ロジスティック回帰で切片項は使用しないこととする。


|sex|height|
|:---|:---|
|男|162|
|女|158|


##問1:ロジスティック回帰を行うためにデータを書き換えよ

ロジスティック回帰では目的の分類クラスを1,0に変換することで、クラスを確率として扱う。
1=100%として考えることで、0~1までの数値の回帰によってクラス分類を達成する。
そのため以下のように書き換えた。
また、切片項を使用しないので、事前に得られている平均値から各データを引いておく。

|sex|height|
|:---|:---|
|1|2|
|0|-2|


##問2：1から0までの値で出力するための変換関数(シグモイド関数)を書け

ロジスティック回帰ではシグモイド関数によって1から0の間に値を変換する。
数式に苦手感を感じる人はexp(x)という記号が2.718のx乗だと考えてくれればいい。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/56bd4797-2f58-4515-fb5b-38f81c9ee454.png)

関数として描画すると以下のような形になる。


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/10c3a5c5-c701-442f-a15a-909bc2e452d6.png)


##問3：ロジスティック回帰の尤度関数から負の対数尤度関数までの流れを記述せよ


まずロジスティック回帰ではシグモイド関数に

- 男の時(sex = 1)のデータ(height = x)を入力すると1に近くなり(大きくなり)
- 男でない時(sex = 0)のデータでは0に近くなる(小さくなる)

という内容を満たすように、係数**w**をデータ**x**にかけることで調節していくのが目的の一つである。

ここで

- **「男の時は、男でないときは」**とかプログラムじゃないんだから一つの式に表して計算したい
- **「大きくする、とか小さくする」**とかどっちかに統一しろや

って気持ちがあるのでこれを満たすような関数を考えてやる。


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/53123809-adf9-ac34-8d9d-66d01f5cd346.png)

そもそも、シグモイド関数にxを入れた時に出てくる値yは0~1の確率値として考えて左辺をP()で表す。
(推定値なので$\hat{y}$や$\hat{sex}$とした方が適切かもしれない)

- 右辺の式はsexが1ならば後ろの部分は0乗となり前部分だけが残る
- sexが0ならば後ろの部分だけが残り、シグモイド関数にxを代入したときに0に近いほど1から引かれる値が大きく成功とみなされる
- この計算を各データに行い、最大になるような**係数w**を求めたい
- 最大を考える時にすべての出力を掛け算して最大になるときのwを考えよう


以上のことからwに関する関数が考えられた。
この関数を**尤度関数**と呼ぶ。

総乗記号を使い一般化して書くと
xに添え字iを書き忘れたが

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/d84adaea-238f-b35a-c6fe-1321f2d2ebd5.png)

という式で表現できる。


しかし、

- 0~1の間の数字を何度も掛け算していくと0に極めて近い値になってしまう
- 掛け算は計算機君でもめんどくさがる

ので、簡単な足し算にしたいよね。
ってことで掛け算の性質を足し算の性質に変換するlog君に登場してもらう。

すると、総乗記号は総和記号に書き換わり以下のように書き直せる。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/149027e6-0ae2-fefe-5dd2-694928c9ea62.png)


尤度関数に対数logを使ったので、これを**対数尤度関数**と呼ぶ

まだまだ変形すんで。
そもそも尤度関数も対数尤度関数も最大にしてくれるような係数wを求めるのが目的だった。
最適化では最小化問題に書き換えるという慣習がある。
~~(きっと計算数値が小さくて済むからだと考えているが違ったらゴメンなさい。教えて。)~~

対数のグラフlog(x)を-log(x)に書き直してみると逆転しているだけだってことがわかると思う。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/c23bde8c-5956-7c29-f6f4-40804b8ce48e.png)


なので全体にマイナスを掛け算してやると


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/fc0a57e6-516a-d31d-05b3-d1a8f379890f.png)

となり、これを**負の対数尤度関数**と呼ぶ。
よく**交差エントロピー**などとも呼ばれる

##問4：負の対数尤度関数を使ってどうやって係数wを決めていくのか

負の対数尤度関数と、

パターン1.男のデータの時の尤度関数の値の関係は

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/9b375f9f-274b-e4a5-4f2b-2b3c5ce1cf9f.png)



パターン2.男でないときの関係

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/22b464a5-d69c-5be6-0644-e366f0c64efd.png)

となる。つまり、負の対数尤度関数の値が小さくなる時は、関数の傾きが小さくなっている時である。
よって関数の傾きを小さく
というか、0にするようにwを決めていけばいいことがわかる。


##問5：じゃあ傾き求めてみて

###5-1全体を微分する準備
傾きを求めるために負の対数尤度関数の導関数を求めていく。

合成関数の微分は複雑なので**みっちり**と順を追って記述していく。

負の対数尤度関数を係数wに関する関数J(w)として考える
あと、sexって書くのが面倒なのでsとして省略した。
微分は各項に対して行うものなので、前の項と後の項を微分する。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/02a1147b-0441-0dad-eda2-19cdaeb92640.png)

###一般規則の復習
ここで対数logの中に関数f(x)が入っている時の微分は、合成関数の基本的な規則より

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/5dcc90ab-9d4b-6c73-5369-95460b7b18b5.png)

対数の微分は逆数となり、括弧内の関数の微分を後から掛け算することで求まると思い出しておく。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/d65b7e07-d6aa-ecd7-e6da-cdefee0d7b75.png)

つまり計算するとこうなる。

###5-2対数の微分後のちょっとした微分処理と式整理

後ろの項の


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/3d6076e5-970f-7e0b-a0ac-394933502063.png)

は計算すると


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/b6b6489f-5d48-cffe-0106-4526b2cf5487.png)

となり、前の項と共通になるのでくくりだせる

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/e05c2ce3-e096-2460-3c07-7f87d77001e6.png)

あと残す微分をやっつけていく。

###5-3シグモイド関数の微分

wxが計算をややこしくしそうなので、zに置き換えてから、合成関数として求めてやる。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/2be012d3-bdeb-8a1f-6dab-719c97e365d6.png)


zの微分は求めるの簡単

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/9ba73872-356f-c23e-fb17-a3666e4b7416.png)

シグモイド関数の微分も合成関数の微分として順番に解いていく

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/2e724f4a-2e37-1e14-08fd-77d64ed60149.png)

なんでこんなこと考えつくんだろって思うけど式変形していくと、

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/ddfa8c92-3d81-1519-7bc9-208a47481231.png)

となる。
最後にはzに関するシグモイド関数が出てきた。

wxをzと置いていたことを思い出して、代入してやると、



![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/03a218b1-25ca-fce8-a6b0-e6a91884c612.png)


zに関するシグモイド関数の微分を計算し終えることができた。

先ほどのZと置いて合成関数の微分を行った結果を合わせると

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/139ff8db-807e-748f-926e-1557a90d6569.png)


となる。

###5-2に5-3の計算結果をあわせる

さて、少し5-2戻って式変形をする
純粋に式を変形しているだけ。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/9ea2ffa2-32f0-a0b9-e41b-2beba43d4b5e.png)

~~どうしてこんな式変形が思いつくんや~~
すると、角括弧内とシグモイド関数の微分によって出てきた部分が約分できる。

あとは純粋に綺麗にしていくと

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/65cf832e-7f8f-094b-4722-c698926125cd.png)

すっごい綺麗にまとまった。

###5-4係数を求め、もともとのロジスティック回帰と、学習後のロジスティック回帰を比較せよ


もともとのシグモイド関数にデータを入力した場合


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/e5fbf738-64a9-f356-6a08-f8597d900430.png)


これでも0.5よりも小さいか大きいかで判別させたら十分機能はしますが、
今回はこれをより0と1に近づけていく学習なので、
勾配降下法によって学習させていきます。

更新式は

>新しい係数 <- 係数 - 学習率*(シグモイド関数(-w*x) - ラベルy)x

です。
勾配降下法については[過去記事](https://qiita.com/Ringa_hyj/items/529ed71407cdec54a160)で紹介しております。

初期値を0として係数を学習させていくと



![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/854ed00c-a61b-14ad-ac01-55dc8ca19d8d.png)



係数=5あたりでゆっくりになってきました。
どれだけ大きくしていってもうまく分けられるのできっと上昇し続けるのではないかと思います。

そこで得られた係数を使ってシグモイド関数を書いてみます。


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/1d8ce86c-0dc3-2041-392f-4f72127fb31d.png)

S字の曲線はより急な変化になり、1や0にうまくデータを振れていることがわかりました。

ということで、係数を求め、シグモイド関数の比較ができました。

ちなみにRのglm関数で求めた場合は係数=11になりました。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/7e756315-7a73-00ac-158e-a28a49aac9a2.png)



#以上
スッキリした。

以下コード

```r
x<-seq(-50,50,0.1)

sigmoid <-function(x){
  A<-1
  B<-1+exp(-x)
  return(A/B)
}

y<-sigmoid(x)
plot(x,y,type="l",main="sigmoid function")

x<-seq(0.1,20,0.1)
y1 <-log(x)
y2<- -log(x)
plot(x,y1,type="l",xlim=c(0,20),ylim=c(-5,5),main="red = -log(x)")
points(x,y2,type="l",col="red")

plot(y,-log(y),type="l",xlab = "sigmoid(-wx)",ylab="cross ent")

x<-seq(-50,50,0.1)
y3<-1-sigmoid(x)
y<-sigmoid(x)

plot(y,-log(y3),type="l",xlab = "sigmoid(-wx)",ylab="cross ent")

#################

x<-c(2,-2)
y<-c(1,0)

plot(x=x,y=y)
data<-data.frame(x=x,y=y)

w_0 <- 0
lr <- 0.5
w_list<-NULL
ite<-10000
lis<-seq(1,ite,length.out = 100)
for(i in 1:ite){
  rand<-sample(2)[1]
  data_x<-data$x[rand]
  label<-data$y[rand]
  w_new <- w_0 - lr * (sigmoid(x=w_0*data_x) - label)*data_x
  w_0 <- w_new
  if(sum(lis%in%i)==1){
    w_list<-c(w_list,w_new)
  }
}
w_0
plot(1:length(w_list),w_list)
#dev.off()

dd<-seq(-3,3,0.1)
plot(dd,sigmoid(dd))
abline(v=data$x[1])
abline(v=data$x[2],col="red")

sigmoid2 <- function(x,w){
  A<-1
  B<-1+exp(-x*w)
  return(A/B)
}

dd<-seq(-3,3,0.1)
plot(dd,sigmoid2(dd,w=w_0),type="l")
abline(v=data$x[1])
abline(v=data$x[2],col="red")

res = glm(y ~ .-1,data, family=binomial)
res_sum<-summary(res)
fit = fitted(res)
plot(data$x, fit,col="red")
points(dd,sigmoid2(dd,w=res_sum$coefficients[1]),type="l")
```