要約
いまさらロジスティック回帰?と思うかもしれないが, もう火鍋の話はしない. 昔書いたやつを読み返したら中途半端だったので改めて(2値)ロジスティック回帰のいくつもある表現について書きたくなった. 昔書いたやつというのは以下のことである.
https://ill-identified.hatenablog.com/entry/2015/03/01/170623
今度は以下の4種類に触れる. 最初の2つは前回の記述を推敲しただけでほとんど同じである.
- 経済学でよく使われるロジットモデルの潜在変数モデルによる表現
- 一般化線形モデル (GLM) の特殊形としてのロジスティック回帰
- 機械学習の分類タスクとしてのロジスティック回帰 (分類) モデル
- 事後確率で見たロジスティック回帰
各教科書に載っている内容を集めただけなので, 目新しい話はない.
はじめに
最もシンプルなモデルであるロジスティック回帰には, いくつもの表現がある. どのモデルでも得られる結果は同一だが, 導出の過程が異なる.
GLM の特殊形としてのロジスティック回帰
久保 (2012) (通称緑本) では, 線形予測子をロジットリンク関数で変換したものとしてロジスティック回帰を解説している. 線形予測子 y_i^\ast を次のように定義する. これ以降と比較しやすいよう, 本紙の記述とは異なる記号を使っている.
y_{i}^{\ast} = \alpha+\beta x_{i}+\gamma z_{i}
x_i, zi はそれぞれ説明変数で, \alpha, \beta, \gamma がパラメータである. GLMでは, 線形予測子を変換する関数をリンク関数はと呼んでおり, ロジットリンク関数とはロジスティックシグモイド関数 (以下, 単にシグモイド関数と呼ぶ)のことである. 線形予測子をシグモイド関数で変換した q_i は以下のようになる.
\begin{align}
q_i &=\mathit{logistic}(y_i^{\ast})\\
&= \frac1{1+exp(−y_i^{\ast})}
\end{align}
となる. シグモイド関数は, (-\infty, \infty) の定義域に対して (0,-1) の値域への変換ができる. つまり, シグモイド関数によって線型予測子を確率として扱える. q_i を y_i=1 の確率とおきかえると, y_i がとる値は二項分布になる. これをもとに尤度関数を作成すると,
\mathcal{L} = \prod_i \binom{N_i}{y_i} q_i^{y_i} (1−q_i)^{N_i−y_i}
となる. しかしながら, 久保本は y_i が 0, 1 以上の値を取る場合の式を書いている. 厳密には, GLMの定義では, 上記の二項分布の式がロジスティック回帰と呼ばれるため, 後述する計量経済学と機械学習での定義とは少し異なる. しかし今回は2値のロジスティック回帰のため, y_i がゼロか1しかない場合だけを考える. N_i も1になるので, 二項分布の最も単純な形式であるベルヌーイ分布となり, 尤度関数はもう少しシンプルに, 以下のように書ける.
\mathcal{L} = \prod_i q_i^{y_i} (1-q_i)^{y_i}
対数尤度関数は以下のように表せる.
\ln \mathcal{L} = \sum_i \left( y_i \ln q_i + (1-y_i) \ln (1-q_i) \right)
潜在変数モデルとしてのロジスティック回帰
計量経済学においては, ロジスティック回帰はロジットモデルと呼ばれることが多い. そして教科書ではロジットモデルを説明する時, 多くの場合で潜在変数 (latent variable) モデルを用いて表現している. 潜在変数とは, 結果を表す変数が直接観察できないという意味である. 経済学では, 人間の意思決定は効用関数で数量的に表現できるが, 多くの場合で効用関数は直接観測できないという想定で理論分析を行っている. 基本的なロジスティック回帰では, 効用は線型回帰と同じように, 説明変数に依存して決定すると仮定している.
y_i^{\ast} =\alpha + \beta x_i + + \gamma z_i
ロジットモデルでは, この効用の大きさが本来の結果変数だが, 効用は観測はできない潜在変数であり, 行動として現れる ゼロか1かの2値の離散変数 y_i を代わりの被説明変数として扱う.
潜在変数である効用を y_i^{\ast} と表すと, 行動として現れる結果変数 y_i のどちらを選ぶかの意思決定が, ゼロを境に決まると仮定する.
y_i= \begin{cases}
1 & \text{if} y_i^{\ast} + \varepsilon_i \geq 0\\
0 & \text{if} y_i^{\ast} + \varepsilon_i < 0
\end{cases}
ここで, \varepsilon_i は誤差項で, 標準ロジスティック分布にしたがうと仮定する. 標準ロジスティック分布とは, 位置パラメータがゼロ, 尺度パラメータが1のロジスティック分布のことで, 累積分布関数は以下で与えられる.
\Lambda(x) = \frac1{1 + \exp(-x)}
よって, y_i=0, 1 である確率, つまり選択確率は, それぞれ以下のようになる.
\begin{align}
\mathrm{P}(y_i = 1) &= \mathrm{P}(y_i^{\ast} + \varepsilon_i \geq 0)\\
&= \mathrm{P} (\varepsilon_i \geq -y_i^{\ast})\\
&= 1 - \frac{1}{1 + \exp(-y_i^{\ast})}\\
&= \Lambda(y_i^{\ast}) \\
\mathrm{P}(y_i = 0) &= \mathrm{P}(y_i^{\ast} + \varepsilon_i < 0)\\
&= \mathrm{P} (\varepsilon_i < -y_i^{\ast})\\
&= \frac{\exp(-y_i^{\ast})}{1 + \exp(-y_i^{\ast})}\\
&= 1 - \Lambda(y_i^{\ast})
\end{align}
潜在変数と説明変数の関係式と合わせると, この選択確率はGLMのロジット q_i と同じであるとわかる.
さらに, ロジットモデルの尤度関数は, 選択確率が背反であるから, 2つの選択確率の積で表せる. よって, 尤度関数もGLMでのロジスティック回帰と同じになる.
\begin{align}
\mathcal{L} &= \prod_i \Lambda(y_i^\ast)^{y_i} \Lambda(y_i^{\ast})^{1-y_i}\\
&= \prod_i q_i^{y_i} (1-q_i)^{1-y_i}
\end{align}
ロジットモデルはこのように, 効用関数の概念を類推しやすい定式化から始まって, ロジスティック回帰と同じ式に帰着している. 効用 y_i^{\ast} が一定以上高ければ, 実際の行動 y_i=1 として現れ, 逆にその行動から得られる効用 y_i^{\ast} が小さければ, y_i=0 という行動, 言い換えるなら「行動をしない」 という選択に現れる. そして, その効用を決定する要素が説明変数である. ある商品を買う (y_i =1) か買わない (y_i=0) かを決めるのは, その商品の価格や品質, あるいは消費者の好みで決まると考えると, これらが説明変数に対応する (これは混合ロジットモデルと呼ばれる). さらに, 潜在変数モデルの考え方は多項選択モデルへ拡張することもできため, 多項ロジスティック回帰でも同じ要領で定式化できる. 詳細は, 昔書いた話 (https://ill-identified.hatenablog.com/entry/2014/06/10/000354) や, より専門的な話が知りたければ, Train (2009), Wooldridge (2010) などがおすすめである.
なお, 個人レベルの意思決定のデータがなく, 大勢の意思決定をいくつかのグループに集計したものだけが使える時に, そのグループ内で何割が選択している, という割合を被説明変数にした回帰モデルを集計ロジットモデルというらしい. 私は知らなかったが, 交通政策の分野なんかで使われることがあるらしい.
機械学習の分類モデルとしてのロジスティック回帰
機械学習でよく用いられる, 分類問題としてマージンから損失関数を構成する定式化を紹介する. 例えば 杉山 (2013), 中川 (2015) では, この定式化で説明されている.
まずは, ロジスティック回帰の定式化の前に, マージンと損失 (loss) の考え方を導入する. ここまでで紹介した2つのロジスティック回帰とは違い, 機械学習の2値分類問題では結果変数を -1, 1 の2通りと定義することが多い. しかし, (0, 1) から (-1, 1) への変換は 2x-1 で簡単にできるため, 本質的な問題とはならない. y_i を正しく分類するモデルを作るためには, 誤まって分類した時は1を, 正しく分類したときは0を返す関数を考える. 最尤法では尤度関数を最大化することでパラメータが決まるが, 機械学習では損失関数を最小化することで決まる.
線形分類モデルでは, マージンは, ラベルと線型回帰モデルの右辺の積で表される. ここで, これまでのロジスティック回帰の代数との混同を避けるため, ラベルを l_i と表すと, マージンは以下のようになる.
\begin{align}
m_i &= l_i \cdot (\alpha + \beta x_i + \gamma z_i)\\
&= l_i y_i^\ast
\end{align}
機械学習では, 分類を誤った時の損失を最小化する問題として定式化する. ラベルは-1か1のどちらかなので, y_i^{\ast} との積であるマージンが正ならば分類が正しいということになるから, 以下のような損失関数が考えられる. これを 0/1損失と呼ぶ.
\begin{align}
\mathit{Loss} &= \sum_i \mathit{loss}(l_i, y_i^\ast)\\
\end{align}
ここで, \mathit{loss} は, 以下で定義される.
\mathit{loss}(l_i, y_i^\ast) =
\begin{cases}
1 & \text{if } m_i > 0\\
0 & \text{if } m_i \leq 0
\end{cases}
しかし, このような損失関数は非連続であり, 組み合わせ最適化問題に属するため, 現実的に計算できない可能性が高い. そこで, 機械学習では0/1損失を近似できて計算しやすい関数で置き換える. そのうち1つがロジスティック損失で近似したロジスティック回帰であり, あるいは2乗損失を利用した最小二乗法であり, ヒンジ損失を使用したサポートベクターマシン (SVM) であり, 指数損失を利用したAdaBoostであったりする. こういった損失関数は代理 (surrogate) 損失という. この中ではロジスティック損失は \mathit{loss} = \ln(1+\exp(m_i)) と定義される.
p(y_i=1\mid x, z) に対応するラベルは l_i=1 なので, 確率は 1/(1+\exp(-y_i^\ast)) = 1/(1 + \exp(-m_i)) と表せる. 同様に y_i=0, つまりl_i=-1 のとき, p(y_i=0\mid x, z) = 1/ (1 + \exp (-m_i)) と表せるので, この表記での対数尤度は以下のように表せる.
\begin{align}
\ln \mathcal{L} &= \ln \prod_i p(y_i=1\mid x_i, z_i)^{y_i}p(y_i=0\mid x_i, z_i)^{1-y_i}\\
&= \ln \prod_i (1+\exp(-m_i))^{-1}\\
&= -\sum_i \ln (1+\exp(-m_i))
\end{align}
対数尤度が, 符号が反転したロジスティック損失となった. 機械学習では損失関数を最小化するから, これは対数尤度の符号を反転して最大化することと同じである. よって, この式から, ロジスティック回帰の対数尤度最大化とロジスティック損失最小化が等価であることがわかる.
事後確率の近似としてのロジスティック回帰
Bishop (2006) の Pattern Recognition and Machine Learning, いわゆる『パターン認識と機械学習』は機械学習の名を冠しているが, 上記のようなマージンを使った説明ではなく, 事後確率を使用して導出している.
\begin{align}
P(y_i=1\mid x) &= \frac{p(x\mid y_i=1) P(y_i=1)}{P(x\mid y_i = 1)p(y_i=1) + p(x\mid y_i=0)p(y_i=0)}\\
&= \frac1{1 + \frac{p(x\mid y_i=0)p(y_i=0)}{p(x\mid y_i=1)p(y_i=1)}}
\end{align}
なお, 上の式では事後確率は本来 x, z で条件付けられるべきだが, 冗長なので z は省略する.
y_i^\ast = \ln \frac{p(x\mid y_1=1)p(y_i=1)}{p(x\mid y_i = 0)p(y_i=0)}
とおくと, 事後確率は以下のようになる.
\begin{align}
p(y_i=1\mid x) &= \frac{1}{1+\exp(-y_i^\ast)}\\
&= \mathit{logistic}(-y_i^\ast)
\end{align}
ここで, p(x\mid y_i) は, データのうち各 y_i に属するものの特徴量 x の条件確率である. この x の条件分布が正規分布だと仮定するなら, 以下で表せる.
p(x\mid y_i=k) = \frac1{\sqrt{2\pi}\sigma_k} \exp\left(-\frac{(x-\mu_k)^2}{\sigma_k^2}\right),\, k = 1, 2
ここでさらに, 各ラベルの標準偏差 \sigma_k が等しく \sigma だと仮定すると, \alpha, \beta は
\begin{align}
\alpha &= -\frac1{2}\frac{\mu_1^2}{\sigma^2} + \frac1{2} \frac{\mu_2^2}{\sigma^2} + \ln \frac{p(y_i=1)}{p(y_i=0)}\\
\beta &= \frac{(\mu_1 - \mu_2)}{\sigma^2}
\end{align}
と表せる. よって, ここでも線形予測子とシグモイド関数の関係が現れる.
\begin{align}
y_i^\ast &= \ln \frac{p(x\mid y_1=1)p(y_i=1)}{p(x\mid y_i = 0)p(y_i=0)}\\
y_i^\ast &= \alpha + \beta x
\end{align}
以降の尤度関数を構成する箇所は, GLMやロジットモデルと同じなので省略する.
この説明は一見すると, 導出仮定がややこしく, 式は複雑で洗練されておらず, いびつに見えるかもしれない. しかし, ここまでの3種類の説明では, 事前確率 p(y_i) にあたる要素に全くふれられていなかった. この説明で初めて事前確率 p(y_i) という概念が登場し, パラメータにどう影響するかが示唆された.
現実では, この事前確率に対処する方法は様々である. 事前確率が各ラベルで均等である, つまり訓練データ内に均等な頻度で登場しているならあまり気にする必要はないが, 不均衡である場合はナイーブベイズモデルであったり, データのリサンプリングによって補正したりといった応用方法がある.
参考文献
Bishop, Christopher M. 2006. Pattern Recognition and Machine Learning. Information Science and Statistics. New York: Springer. https://www.microsoft.com/en-us/research/people/cmbishop/prml-book/.
Train, Keneth. E. 2009. Discrete Choice Methods with Simulation. 2nd ed. New York: Cambridge University Press. https://eml.berkeley.edu/books/choice2.html.
Wooldridge, Jeffrey M. 2010. Econometric Analysis of Cross Section and Panel Data. 2nd ed. The MIT Press. https://mitpress.mit.edu/books/econometric-analysis-cross-section-and-panel-data.
中川裕志. 2015. 機械学習. Edited by 東京大学工学教程編纂委員会. 東京大学工学教程 情報工学. 丸善出版.
久保拓哉. 2012. データ解析のための統計モデリング入門. 岩波書店.
杉山将. 2013. イラストで学ぶ機械学習: 最小二乗法による識別モデル学習を中心に. 東京: 講談社.
Discussion