背景
「ディープラーニングと物理学」という本面白いです。
分類問題を統計力学の問題として設定するとシグモイド（としか書いてないけどロジスティック回帰）が出てくるあたりが特に。
ただ、カノニカル分布ってどんな仮定を使ってるんだっけとか思い出せず。。。
最尤推定とかベイズ統計とか正直ピンときてないので、とても簡単な問題でいろいろ考えてみました。
問題設定
外場（磁界とか）の入力$\vec{x}$ があったときに対象の状態（スピンとか）の観測結果$d$が0と１のどちらになるかを当てるためのモデルを作りたい。
実験はたくさん（N回）していてデータがある。
i番目の実験結果はこのようになっている。
入力 : $\vec{x}[i] \in \bf{R}^2$
出力： $d[i] \in \lbrace 0, 1 \rbrace $
そして、線形分離できる感じがする実験結果があったとする。
グラフの横軸が$\vec{x}$のx成分で縦軸がy成分。
ドットの色が$d[i]$の値(0 or 1)を表す。
まずは普通にロジスティック回帰
式の展開は「はじパタ」の6.4.3章で入力をベクトルにしているだけです。
なんでこんなことをしているかの「お気持ち」を可能な限り書き出してみています。
「2クラスロジスティック回帰モデルのパラメータ$\vec{w}$と$b$の最尤推定を考える。」
問題設定とすぐ上の「」内を合わせて言い換えると、
-
今回のＮ回の実験で結果（$\vec{x}[i]$のときに出力$d[i]$になる）が出た。
-
モデルはロジスティック回帰を仮定しよう。
（最適だって根拠はないけど線形分離できそうな実験結果だし。）
$d[i]$ が1になる確率が$\pi_i$、0になる確率を$１-\pi_i$として$\pi_i$は以下だとするということ。$$\pi_i = \sigma (\vec{w} \cdot \vec{x}+b) = \frac{e ^{\vec{w} \cdot \vec{x}+b}}{1+e ^{\vec{w} \cdot \vec{x}+b}}$$
-
モデルのパラメータ $\vec{w}$ と $b$ を決めたい。
（けど、ここまでの事実と仮定だけじゃ決められない。不良問題だ。）
今回のN回分の実験結果が出る確率を最大にするようなパラメータが設定されていたのだろう。
（必然性はないけど、この仮定を使うことにする。）
と考えて、計算を始めます。
Nの試行の尤度関数（N回のの試行が実験結果と同じになる確率）は、この式。
$$L(\pi_1,...,\pi_N)=\prod_{i=1}^{N}f(d[i]|\pi_i)=\prod_{i=1}^{N}{\pi_i^{d[i]}(1-\pi_i)^{(1-d[i])}}$$
ここで$f(d[i]|\pi_i)$は、$d[i]$が1のとき${\pi_i}$になって、$d[i]$が0のとき$(1-\pi_i)$になる。よくやるやつ。
尤度関数$L$最大化する代わりに負の対数尤度関数を最小にする問題にする。交差エントロピーの形になる。
$$\mathcal{L} = -\ln L = -\sum_{i=1}^{N}\left({d[i]\ln\pi_i + (1-d[i])\ln(1-\pi_i)}\right)$$
$\pi_i$を代入して整理する。
$$\mathcal{L} = -\sum_{i=1}^{N}\left(d[i](\vec{w} \cdot \vec{x}+b) -\ln(1+e ^{\vec{w} \cdot \vec{x}+b})\right)$$
最小値を求めたいので偏微分した式が0になるのが解。
すごくきれいな式だけど解析的には解けなくて数値計算になる。
$$\frac{\partial\mathcal{L}}{\partial\vec{w}}= -\sum_{i=1}^{N}\left({d[i]\vec{x}- \frac{\vec{x}e ^{\vec{w} \cdot \vec{x}+b}}{1+e ^{\vec{w} \cdot \vec{x}+b}}}\right) = \sum_{i=1}^{N}\vec{x}\left(\pi_i-d[i]\right)$$
$$\frac{\partial\mathcal{L}}{\partial b}= -\sum_{i=1}^{N}\left({d[i]- \frac{e ^{\vec{w} \cdot \vec{x}+b}}{1+e ^{\vec{w} \cdot \vec{x}+b}}}\right) = \sum_{i=1}^{N}\left(\pi_i-d[i]\right)$$
数値計算してみる。
統計力学になかなかたどり着きませんが、数値計算してみます。
ここまで式を出しておきながら、計算はscikit learn で、
clf = LogisticRegression(random_state=0).fit(X, y)
clf.predict(X[:2, :])
では寂しいものがあります。
L-BFGS法を直接使って計算します。
この方法で先ほど解析的に出した$\frac{\partial\mathcal{L}}{\partial\vec{w}}$や$\frac{\partial\mathcal{L}}{\partial b}$が役立ちます。
データはこちらから頂きました。
http://arduinopid.web.fc2.com/N55.html
L-BFGS法のライブラリとしてはliblbfsを使わせてもらいました。
詳細省きますが、C言語の実装をこちらに置いています。
82～84行目が$\frac{\partial\mathcal{L}}{\partial\vec{w}}$と$\frac{\partial\mathcal{L}}{\partial b}$です。
こちらが計算結果です。dが教師データ。学習後のpが出力です。
データを正規化しない雑なやり方でしたが、3つのパラメータ$w_x$, $w_y$, $b$だけで学習用データをうまく分離できたようです。
Let's check the results!
x = 41, y = 104, d = 0, p = 0.000000 OK
x = 76, y = 244, d = 0, p = 0.000000 OK
x = 168, y = 224, d = 0, p = 0.000000 OK
x = 235, y = 139, d = 0, p = 0.000000 OK
x = 133, y = 115, d = 0, p = 0.000000 OK
x = 489, y = 200, d = 1, p = 1.000000 OK
x = 435, y = 40, d = 1, p = 1.000000 OK
x = 520, y = 213, d = 1, p = 1.000000 OK
x = 248, y = 223, d = 0, p = 0.000000 OK
x = 411, y = 168, d = 1, p = 1.000000 OK
x = 499, y = 182, d = 1, p = 1.000000 OK
x = 444, y = 137, d = 1, p = 1.000000 OK
x = 301, y = 135, d = 1, p = 0.999998 OK
x = 414, y = 200, d = 0, p = 0.000002 OK
x = 418, y = 158, d = 1, p = 1.000000 OK
x = 150, y = 25, d = 1, p = 1.000000 OK
x = 512, y = 194, d = 1, p = 1.000000 OK
x = 394, y = 216, d = 0, p = 0.000000 OK
x = 91, y = 89, d = 0, p = 0.000000 OK
x = 278, y = 64, d = 1, p = 1.000000 OK
統計力学の場合
こちらの教科書を参考にしています。
まずエントロピーと温度を定義します。
その後カノニカル分布を考えます。
物理量であることをほとんど気にしなくてよい論理展開になっているのがポイント。
エントロピーの定義
AとBからなる合計エネルギーがEの孤立したシステムを考える。
Aのエネルギーが$E_A$であるときの状態の数を$Ｗ_A(E_A)$と書くとする。
例えば、同じ速さの粒子でも方向が違えば違う状態なので、それだけ状態の数がある。
ここで、等確率の原理を仮定する。統計力学の超重要な仮定。
「孤立したマクロな物体では、十分に長い時間で見ると、実現可能な量子状態はすべて等しい確率で実現する。」
よって、あるエネルギー配分の実現する確率$P(E_A, E_B)$は以下。
$$ P(E_A, E_B) = \frac{W_A(E_A)W_B(E_B)}{\sum_{E'_A+E'_B=E}W_A(E'_A)W_B(E'_B)} $$
エントロピーは以下のように定義する。一番右がいかにも情報量の形。
$k_B$はボルツマン係数。
$$S_A(E_A)=k_B \log W_A(E_A)=-k_B \log \frac{1}{W_A(E_A)}$$
Aのエネルギーが$E_A$でBのエネルギーが$E_B$である時、状態の数$W_A(E_A)W_B(E_B) $なので、
そのときのエントロピー$S(E_A, E_B)$は以下のように書ける。
$$S(E_A, E_B) = k_B \log W_A(E_A)W_B(E_B) ＝ S_A(E_A) + S_B(E_B) $$
温度の定義
以下を温度の定義にする。
$$\frac{dS(E)}{dE}=\frac{1}{T}$$
これだけだと納得感が無いので、脱線して物理的な意味の話につなげる。
二つの系がつながったAとBが平衡状態となって観測できるのは、確率$ P(E_A, E_B)$が最大になるような$E_A$, $E_B$の組み合わせのときだと考える。
つまりパラメータ$E_A$, $E_B$の最尤推定をする。
$ P(E_A, E_B)$の分母は和をとったあと定数になるので、$ S(E_A, E_B)$ が最大になるときに$P(E_A, E_B)$が最大になる。
エネルギー一定であることに気を付けて、$ S(E_A, E_B)$を微分した結果が0のときが最大値になる。
$$\frac{dS(E_A, E_B)}{dE_A} ＝ \frac{dS_A(E_A)}{dE_A} - \frac{dS_B(E_B)}{dE_B} =0$$
もしも、$\frac{dS_A(E_A)}{dE_A} $が$\frac{dS_B(E_B)}{dE_B} $より大きいと、$E_A$による微分が正なので、$E_A$が増えて$E_B$が減ることで平衡状態に行く。
もしも、$\frac{dS_A(E_A)}{dE_A} $が$\frac{dS_B(E_B)}{dE_B} $より小さいと、$E_A$による微分が負なので、$E_A$が減って$E_B$が増えることで平衡状態に行く。
エントロピーをエネルギーで微分したものの大小関係でエネルギーが流れて平行状態に向かっている。これは温度に似ている。
ということで、以下が温度の定義だとしているのが妥当だと思われる。
$$\frac{dS(E)}{dE}=\frac{1}{T} $$
こうするとT_A < T_BのときにBからAにエネルギーが流れることになる。
今の議論だけだと$\frac{dS(E)}{dE}$はTの単調減少の関数であればよい。
みんなが絶対温度として使ってみて便利だったのはこれ$\frac{dS(E)}{dE}=\frac{1}{T} $だったらしい。
そもそも、熱平衡の状態を推定するのに最尤推定使うのは妥当か？というのは仮定次第。
$W(E)$が$E$により指数関数的に増加する関数であることを仮定できると、平衡状態周辺にしか分布が無い状況になり妥当。
カノニカル分布
注目している系Aのある量子状態のエネルギーがE_nとする。系Bは系Aと比べてとても大きい熱浴と呼ばれる状態だとする。
系Aがエネルギー$E_n$のある状態nにある確率は、対応するBの状態の数$W_B(E_T-E_n)$に比例する。
$$P_n \propto W_B(E_T-E_n) $$
エントロピーの定義$S_B(E)=k_B \log W_B(E)$を使うと、以下のようになる。
$$P_n \propto \exp[\frac{1}{k_B}S_B(E_T-E_n)] $$
ここで$E_n \ll E_T$を使ってテーラー展開し、
$$S_B(E_T-E_n) \approx S_B(E_T) - \frac{dS_B}{dE}_{E= E_T}E_n $$
代入。温度Tの定義$\frac{dS(E)}{dE}=\frac{1}{T} $も使う。
$$P_n \propto \exp\left
[\frac{1}{k_B}\left\lbrace S_B(E_T)-\frac{E_n}{T}\right\rbrace \right] $$
$n$にかかわる項だけ残して、確率を規格化（$\sum_nP_n=1$）するとこの形になる。
$$P_n = \frac{e^{-E_n/k_BT}}{\sum_{n'} e^{-E_{n'}/k_BT}} $$
分類問題のハミルトニアン
こちらの本に戻ってきました。
エネルギーの関数であるハミルトニアンを仮定します。
実験結果を見ると線形識別できそうだし、$d$の一次の項の単純なものにします。
$$ H_{\vec{w},\vec{x},b}(d) = -(\vec{w} \cdot \vec{x}+b)d$$
外場$\vec{x}$のときの条件付確率分布を求める。
カノニカル分布で出した式のエネルギーにハミルトニアンを代入する。
$d$が0か1しか取れなかったことを思い出すと条件付確率は、以下。
$$Q_{\vec{w},b}(d|\vec{x}) = \frac{e^{-H(d)/k_BT}}{\sum_{d=0}^1 e^{-H(d)/k_BT}} = \frac{e^{(\vec{w} \cdot \vec{x}+b)/k_BT}}{1 + e^{(\vec{w} \cdot \vec{x}+b)/k_BT}}$$
真ん中の式が、温度付きソフトマックスの形になっているのも注目ポイント。ちゃんとほしい位置に温度があります。
最後に$k_BT$を$\vec{w},b$に吸収させると最初のロジスティック回帰と同じ形になります。
$$Q_{\vec{w},b}(d|\vec{x}) = \frac{e^{\vec{w} \cdot \vec{x}+b}}{1 + e^{\vec{w} \cdot \vec{x}+b}} = \sigma (\vec{w} \cdot \vec{x}+b)$$
KLダイバージェンス
ハミルトニアンを使った方法では最尤推定を使っていませんでした。
ベイズ統計の立場をとっているから？KLダイバージェンスの最小化を目指しています。
$$Q_{\vec{w},b}(d|\vec{x}) = \sum_{\vec{x},d}P(d,\vec{x})\log\frac{P(d,\vec{x})}{Q_{\vec{w},b}(d,\vec{x})}
= \sum_{\vec{x},d}P(d,\vec{x})\log\frac{P(d|\vec{x})P(\vec{x})}{Q_{\vec{w},b}(d|\vec{x})P(\vec{x})}$$
$$Q_{\vec{w},b}(d|\vec{x}) = -\sum_{\vec{x},d}P(d,\vec{x})\log Q_{\vec{w},b}(d,\vec{x}) + {（\vec{w},b に依存しない項）}$$
データによる経験確率で近似。
$$Q_{\vec{w},b}(d|\vec{x}) \sim = - \sum_{i =1}^N \frac{1}{N}\log Q_{\vec{w},b}(d[i]|\vec{x}[i]) = - \frac{1}{N}\left({d[i]\ln\pi_i + (1-d[i])\ln(1-\pi_i)}\right) = \mathcal{L} / N $$
結局、同じ交差エントロピーにたどり着きます。
結局何を仮定していたか？
統計力学バージョンでは何を仮定していたかというと、、、
- 等重率の原理から確率分布： $ P(E_A, E_B) = \frac{W_A(E_A)W_B(E_A)}{\sum_{E'_A+E'_B=E}W_A(E'_A)W_B(E'_A)} $
- エネルギー保存 ： $E_B = E_T - E_A$
- 熱浴の存在 ：$E_A << E_B$
- ハミルトニアンが線形（モデルの形）
でした。
温度のところで最尤推定しましたが、これは温度の定義の納得感のためなのでノーカウントにしました。
思ったより少ない仮定で導かれていますが、「等確率の原理」はなかなか大きな仮定です。
等重率の原理はベイズ統計でいうところの「事前分布」にあたります。
ベイズ統計の「事前分布」といえば、「人によって設定が変わる困ったやつ」といった話が出てきますが、物理の場合は実験との一致でコンセンサスを取っているようです。
例えば、「黒体輻射」、「気体の比熱」というのは量子力学によって等重率の適用方法が改善されました。
ベイズ統計と統計力学
さらに「ベイズ統計と統計力学は数理的に等価である．」という話があるようです。
引用されているこの本では、分配関数とか自由エネルギーとか物理っぽい用語が飛び交うが、統計力学とは言ってない。
物理よりの説明。
ベイズ予測分布 事後確率最大化推測 最尤推測
参考