車輪の再発明としての二項ロジスティック回帰 feat. 指数分布族

Motivation

ちょっと複雑なモデルを開発・実装しなくてはならないので、単体テスト的な意味で車輪の再発明をする。
サービス提供の質や速度の観点からすると車輪の再発明は望ましくないが、学習者は車輪の再発明をすることを恐れてはいけない、ってどっかの偉人が言っていた気がする（言っていないかもしれない）。

指数分布族に属する分布は対数尤度関数の偏微分（＝スコア関数）が、解析的に求められることができるし、それを利用してGLM（といえる統計モデル）は高速に最適化できるよって話。
新たなモデルを思いついた時、最尤関数だけを求めるのは比較的簡単なのだが、実用性（計算スピード）を考えると勾配も求められないと厳しいのである。

参考資料

大阪大学の舟木剛教授の授業資料（応用システム工学）
高井・星野・野間（2016）『欠測データの統計科学』
WEBに落ちてたプリンストン大学の講義資料 by Barbara Engelhardt（機械学習とかやってる人っぽい）

欠測データの統計科学――医学と社会科学への応用 (調査観察データ解析の実際第1巻)

作者: 星野崇宏,岡田謙介
出版社/メーカー: 岩波書店
発売日: 2016/04/20
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

指数分布族とは

$f( y| \theta ) = exp(a(y)^T b(\theta) + c(\theta) + d(y))$

という形でその密度関数を表現可能な確率分布は、指数分布族（exponential family）である。
ここで $a(y)$ は、確率変数のみで構成された関数のベクトルであり、 $b(\theta)$ は、パラメータのみで構成された関数のベクトル、 $c(\theta)$ はパラメータのみで構成された関数、 $d(y)$ は確率変数のみで構成された関数である。

正規分布をはじめ、統計学で頻出する有名な確率分布はその多くが指数分布族である。

指数分布族の便利な性質

指数分布族にはいくつか便利な性質があるが、その最たるものは対数尤度関数の微分 *1を解析的に求められることである。

対数尤度関数を $l(\theta)$ , その（偏）微分を $l'(\theta)$ と表すこととすると、

$l'(\theta) = a(y)^T b'(\theta) + c'(\theta)$

と書ける。
よく知ってる形の密度関数をわざわざまどろっこしい形に書き直すことの利点はここにある。
比較的用意に勾配の関数が特定できるので、（最尤推定等の）最適化計算を準ニュートン法などで高速に行える。

他にもいくつか便利な性質はあるが割愛（詳しくは舟木先生のスライドにて）。

指数分布族の性質とGLM

さて、我々が慣れ親しんでいるgeneralized linear model（GLM, 一般化線形モデル）は、

従属変数が従う確率分布 $f(y | \theta )$
独立変数の影響をとらえる線形予測子 $\mu = \beta X$
線形予測子と確率分布上のパラメータとの関係を表すリンク関数 $\eta(\mu) = \theta$

といった三つの部分からなる。
我々が知りたいのは、 $\beta_j$ の最尤推定値なので、用いたい勾配というのは $\frac{\partial l}{\partial \beta_j}$ である。

ここで、
$\frac{\partial l}{\partial \beta_j} = \frac{\partial l}{\partial \theta} \frac{\partial \theta}{\partial \mu} \frac{\partial \mu}{\partial \beta_j}$

となるが、前述の指数分布族の性質を利用すれば第一項は解析的に容易に求められる。
また第三項は言うまでもなく $\frac{\partial}{\partial \beta_j}\Sigma_{i=1}^{K} \beta_i x_i = x_j$
第二項は（逆）リンク関数を偏微分すればよい。

かくして、GLMの勾配は容易に求められる。

指数分布としてのベルヌーイ分布

道具立ては整ったので、実装に向けて具体的に二項logistic回帰モデルについてかんがえていく。
上で書いたGLMの三要素に沿って書くと、

従属変数はベルヌーイ分布にしたがう
線形予測子とベルヌーイ分布のパラメータは、logit変換によってリンクされる

というのがlogit modelである。

そこでまず、ベルヌーイ分布を正規分布族として表してみると、
$f(y) = \mu^y (1-\mu)^{1-y} = exp{(y log{\mu} + (1-y)log{(1-\mu)} )= exp(y log{\mu}-ylog(1-\mu) + log(1-\mu))}$
となるので
$a(y) = \begin{pmatrix} y & y \end{pmatrix}^T$
$b(\theta) = \begin{pmatrix} log{\mu} & -log(1-\mu) \end{pmatrix}^T$
$b'(\theta) = \begin{pmatrix} \frac{1}{\mu} & \frac{1}{1-\mu} \end{pmatrix}^T$
$c(\theta) = log(1-\mu)$
$c'(\theta) = -\frac{1}{1-\mu}$

したがって、上の公式を利用して、
$l'(\theta) = a(y)b'(\theta) + c'(\theta) = y\frac{1}{\mu(1-\mu)} -\frac{1}{1-\mu}$
となる。ここから、
$\frac{\partial l}{\partial \beta_j} = \frac{\partial l}{\partial \theta} \frac{\partial \theta}{\partial \mu} \frac{\partial \mu}{\partial \beta_j}= (y\frac{1}{\mu(1-\mu)} -\frac{1}{1-\mu})\mu(1-\mu)x_j$
$= (y_i - \mu_i)x_{ij}$
が導かれる。これを全ケースに関して足し合わせたものを勾配とすればよい。

Rによる実装

以上勾配を求めることができたので、これで車輪の再発明をしてみる。

まず、勾配を返す関数を組む（第一引数がパラメータ、第二引数が従属変数、第三引数が独立変数）

#############------ 勾配関数 ------#############
logreg_grad <- function(prms , y , x, decreasing=F){

#- - 線形予測子の列をつくる- - #
x_mat <- cbind(1,  as.matrix(as.data.frame(x))) #切片項目のため
y <- as.numeric(y)

# N × 1の線形予測子列をつくる
eta <-  as.numeric(matrix(prms , nrow=1 )) %*% t(x_mat)
mu <- 1 /  ( 1+ exp(-1*eta) ) #逆logit 変換

#--- 勾配を計算して返す- - #
y_mu <-  matrix( rep( (y-mu) , ncol(x_mat) ) ,ncol=ncol(x_mat) ) #

grad_mat <- y_mu * x_mat
grad_vec <- apply( grad_mat , 2, sum) 

# y-mu: N × 1 , x_mat = N * K 
return(grad_vec)
} #func

つぎに、尤度関数を組む

#############------ 尤度関数 ------#############
logreg_llik <- function(prms , y , x, decreasing=F){

#- - 線形予測子の列をつくる- - #
x_mat <- cbind(1,  as.matrix(as.data.frame(x))) #切片項目のため
y <- as.numeric(y)

# N × 1の線形予測子列をつくる
eta <-  as.numeric(matrix(prms , nrow=1 )) %*% t(x_mat)
mu <- 1 /  ( 1+ exp(-1*eta) ) #逆logit 変換

# 対数尤度のデータ列をつくる
llik <-  sum( y * log(mu) + (1-y)*log(1-mu))
if(decreasing){
llik <- llik * -1 
} #if

return(llik)
} #logreg_llik(func)

これらを使ってlogistic regression（最尤推定）する関数を組む

############ 上の尤度関数、勾配関数をつかって、logistic regressionする関数 ############
logreg2 <- function(y_var ,x_var, tdf){

y2 <- tdf[[y_var]]
x2 <- tdf[, x_var]
tdf2 <-  as.data.frame(na.omit( cbind( y2, x2)) )
y <- tdf2[,1]
x <- tdf2[ , (2:ncol(tdf2))]

#- - -初期値の設定- - -#
prm0 <- rep(0 , ncol(x)+1 ) 

opt <- optim(par = prm0 , fn = logreg_llik ,gr=logreg_grad , y=y , x = x, control=list(fnscale = -1, maxit=1000,abstol=1e-10),method="BFGS")
return(opt)
} #function logreg2 . . .

#- - 実行- - -#

res1 <- logreg("y" , c("x1","x2"),test_df) #勾配を使わないver

#  計算速度
# ユーザ   システム       経過  
#     0.110      0.004      0.113 

res2 <- logreg2("y" , c("x1","x2"),test_df) #勾配を使うver

#  計算速度
#ユーザ   システム       経過  
#     0.049      0.004      0.052