2019-01-07

Kullback-Leibler Divergenceについてまとめる

統計学機械学習 Python

1. KL-divergenceとは?
2. 諸量との関係
3. 参考書籍

!! お知らせ（2020.06.10） * こちらの記事の英語版を公開しました．よければご覧ください．
Here is the english translation of this post. Please check it if you want.

1. KL-divergenceとは?

　統計学や情報理論をはじめとした広い分野で、KL-divergence*1はたびたび登場します。KL-divergenceは、「尤度比（尤もらしさを比較する尺度）を log 変換し（乗算操作を線形結合に直す、凸関数だから最適化との相性良い）、期待値（確率密度の重み付きの積分、ルベーグ積分）をとったもの」として定義されます。

1.1 定義

　一般に、確率分布 $P, ~ Q$ が確率密度関数 $p(x), ~ q(x)$ をもつとき、KL-divergence（Kullback-Leibler divergence）は、以下のように定義されます。

\begin{eqnarray} D_{KL}( Q \mid\mid P ) := \int q(x) \log \frac{q(x)}{p(x)} ~dx \tag{1} \end{eqnarray}

1.2 基本的な性質

　KL-divergenceは以下のような性質をもちます。

（非負性）非負の値域をもつ。

\begin{align} 0 \leq D_{KL}( Q \mid\mid P ) \leq \infty \tag{2.1} \end{align}

（完備性）値が 0 のとき、2つの確率分布 P と Q は同等。

\begin{align} D_{KL}( Q \mid\mid P ) = 0 ~~ \Leftrightarrow ~~ P = Q \tag{2.2} \end{align}

（非対称性）値は P と Q に対して対称性を持たない。

\begin{align} D_{KL}( Q \mid\mid P ) \neq D_{KL}( P \mid\mid Q ) \tag{2.3} \end{align}

（絶対連続性）値が発散しない限り、Q は P に関して絶対連続である。

\begin{align} D_{KL}( Q \mid\mid P ) \lt \infty ~~ \Rightarrow ~~ P \gg Q \tag{2.4} \end{align}

　例として、2つの正規分布間でKL-divergenceを計算*2すると、以下のような値になります。分布のズレが大きくなるほど、値が増加していることがわかります。

f:id:yumaloop:20190109162028p:plain:w500

1.3 KL-divergenceは距離なのか？

　KL-divergenceは、確率や情報を考える際にとても重要な量であるため、分野や文脈によって様々な名称で呼ばれます。

"KL-divergence"
"KL-metrics"
"KL-information"
"Information divergence"
"Information gain"
"Relative entropy"

　KL-divergenceはつねに非負の値をとるため、これは確率分布 P と Q が存在する空間における距離(metrics)を示していると解釈することができます。しかし、KL-divergenceは次の距離の公理のうち「非負性」と「完備性」しか満たさないため、厳密には距離(metrics) ではありません。

距離 $d(~)$ の公理

非負性　　　　 $d(x, ~ y) \geq 0$

完備性　　　　 $d(x, ~ y) = 0 ~~ \Leftrightarrow ~~ x = y$

対称性　　　　 $d(x, ~ y) = d(y, ~ x)$

三角不等式　　 $d(x, ~ y) + d(y, ~ z) \geq d(x, ~ z)$

　例えば、ユークリッド距離、二乗距離、マハラビノス距離、ハミング距離などは、この4つの条件をすべて満たすため、明らかに距離(metrics)として扱うことができます。一方で、KL-divergenceは距離(metrics)ではなくdivergenceです。"divergence"とは、距離の公理における４つの条件のうち「非負性」「対称性」のみを採用したものであり、"距離(metrics)"の拡張概念です。距離の公理による制約を減らし、抽象度の高い議論をしようというわけです。

　なお、このdivergenceという語は発散と訳されることが一般的で、例えば物理学ではベクトル作用素 div として登場します。divergenceのイメージに該当する日本語はありませんが、相違度、分離度、逸脱度、乖離度などが使われているようです。

　実例として、2つの正規分布 $N(0, 1)$ （青）と $N(1, 2)$ （赤）の間のKL-divergenceを測ってみましょう。左が青からみた赤とのKL距離、右が赤からみた青とのKL距離となっており値が異なることがわかります。なお「等分散正規分布間のKL-divergence」という特殊な条件下では、対称性が成り立ちます。

また，確率分布間の近さを測る尺度としては，KL情報量の他にも次のようなものが知られています．

確率分布 $q(x), p(x)$ 間の"距離"を測る量

${\chi}^2(q ; p) := \sum_{i=1}^{k} {\large \frac{{(p_i - q_i)}^2}{p_i} }$ 　　　　　　　　 ${\chi}^2$ 統計量

$L_1(q ; p) := \int | q(x) - p(x) | ~ dx$ 　　　　　　 $L_1$ -ノルム

$L_2(q ; p) := \int { \{ q(x) - p(x) \} }^{2} ~ dx$ 　　　　　 $L_2$ -ノルム

$I_K(q ; p) := \int { \{ \sqrt{ q(x) } - \sqrt{ p(x) } \} }^{2} ~ dx$ 　　ヘリンジャー距離

$\mathbb{D}(q ; p) := \int f \left( {\large \frac{q(x)}{p(x)} } \right) q(x) ~ dx$ 　　　　　　 $f$ -ダイバージェンス　

$I_{\lambda}(q ; p) := \int \left\{ { \left( {\large \frac{q(x)}{p(x)} } \right) }^{\lambda} - 1 \right\} q(x) ~ dx$ 　　一般化情報量(Kawada, 1987)

$\mathbb{D}_{KL}(q ; p) := \int \log \left( {\large \frac{q(x)}{p(x)} } \right) q(x) ~ dx$ 　　　　 KL情報量

2. 諸量との関係

2.1 KL-divergenceと相互情報量

　エントロピー $H(X)$ 、結合エントロピー $H(X,Y)$ 、条件つきエントロピー $H(X|Y)$ 、相互情報量 $MI(X,Y)$ は、確率密度 $Pr(~)$ を用いて以下のように定義されます。*3

\begin{align} H(X) & := - \int Pr(x) \log Pr(x) ~dx \tag{3.1} \\ H(X,Y) & := - \int Pr(x,y) \log Pr(x,y) ~dy~dx \tag{3.2} \\ H(X|Y) & := - \int Pr(x,y) \log Pr(x|y) ~dx~dy \tag{3.3} \\ MI(X,Y) & := \int \int Pr(x,y) \log \frac{Pr(x,y)}{Pr(x)Pr(y)} ~dxdy \tag{3.4} \end{align}

2つの確率変数 X, Y について、相互情報量 MI（: Mutual Information）によって相互の関係性が明示されます。

f:id:yumaloop:20190107163949p:plain:w400

\begin{align} MI(X,Y) & = H(X) - H(X|Y) \tag{3.5.1} \\ & = H(Y) - H(Y|X) \tag{3.5.2} \\ & = H(X) + H(Y) - H(X,Y) \tag{3.5.3} \end{align}

　ここで、KL-divergenceと相互情報量 MI には以下の関係が成り立ちます。つまり、相互情報量 MI(X, Y) は、「確率変数 X と Y が独立でない場合の同時分布 P(X, Y)」と「独立である場合の同時分布 P(X)P(Y)」の離れ具合（平均的な乖離度合い）を示しているのだと解釈できます。

\begin{align} MI(X, Y) & = D_{KL} \bigl( Pr(x, y) \mid\mid Pr(x)Pr(y) \bigr) \tag{3.6.1} \\ & = \mathbb{E}_{Y} \bigl[ D_{KL} \bigl( Pr(x|y) \mid\mid Pr(x) \bigr) \bigr] \tag{3.6.2} \\ & = \mathbb{E}_{X} \bigl[ D_{KL} \bigl( Pr(y|x) \mid\mid Pr(y) \bigr) \bigr] \tag{3.6.3} \\ \end{align}

（例）3.6.2の式展開

\begin{eqnarray} MI(X,Y) & = & \int \int Pr(x,y) \log \frac{Pr(x,y)}{Pr(x)Pr(y)} ~dxdy \\ & = & \int \int Pr(x|y)Pr(y) \log \frac{Pr(x|y)Pr(y)}{Pr(x)Pr(y)} ~dxdy \\ & = & \int Pr(y) \int Pr(x|y) \log \frac{Pr(x|y)}{Pr(x)} ~dx~dy \\ & = & \int Pr(y) \cdot D_{KL} \bigl( Pr(x|y) \mid\mid Pr(x) \bigr) ~dy \\ & = & \mathbb{E}_{Y} \bigl[ D_{KL} \bigl( Pr(x|y) \mid\mid Pr(x) \bigr) \bigr] \end{eqnarray}

2.2 KL-divergenceと対数尤度比

　ベイズ推定や統計モデリングでは、しばしば「真の分布 $q$ を、 $p(x|\hat{w})$ （パラメータの推定値 $\hat{w}$ と確率モデル $p(x|w)$ の組み合わせ）推定する」というモチベーションが起こります。そのため、２つの分布のズレを測りたいとき、あるいは推定誤差を損失関数やリスク関数に組み込んでパラメータに対する最適化問題をときたいときに、KL-divergenceを使います。

　KL-divergenceは最尤法をベースに、尤度比検定（Likelihood ratio test）、ベイズ因子（Bayse factor）、赤池情報量基準（AIC, Akaike's Information Critation）などのモデル選択手法*4と深い関わりを持ちます。

　真の分布 $P$ に対する推定分布 $p$ のKL-divergence： $D_{KL}(Q \mid\mid P)$ は次式のように、「 $q$ と $p$ の対数尤度比に対して、真の分布 $q$ で平均操作をおこなった値」です。

\begin{align} \left( 対数尤度比 \right) = \log \frac{q(x)}{p(x)} \tag{4} \end{align}

\begin{eqnarray} D_{KL}( Q \mid\mid P ) & := & \int q(x) \log \frac{q(x)}{p(x)} ~dx \\ & = & \mathbb{E}_{X} \left[ \log \frac{q(x)}{p(x)} \right] \left(: 平均対数尤度比 \right) \tag{5.1} \\ \end{eqnarray}

　　モデル比較・選択の評価値としてKL-divergenceを用いる場合、次式のように「KL-divergenceを最小にすること」と「平均対数尤度を最大にする」ことは等価になります。

\begin{eqnarray} D_{KL}( Q \mid\mid P ) & = & \mathbb{E}_{X} \bigl[ \log q(x) \bigr] - \mathbb{E}_{X} \bigl[ \log p(x) \bigr] \tag{5.2} \\ & \propto & - \mathbb{E}_{X} \bigl[ \log p(x) \bigr] \left(: - 平均対数尤度 \right) \tag{5.3} \end{eqnarray}

パラメトリックな確率モデル $f(x|\theta)$ （例えば線形回帰モデル）を考えた場合、何らかの損失関数 $L$ に対して最適なパラメータ $\theta_{0}$ と、その推定値 $\hat{\theta}$ に対して、確率モデルの誤差は以下のようにKL-divergenceで表現されます。（ただし、 $\ell( \cdot| x )$ を対数尤度関数とする。）

\begin{align} \left( 対数尤度比 \right) = \log \frac{f(x|\theta_{0})}{f(x|\hat{\theta})} \tag{6} \end{align}

\begin{eqnarray} \hat{\theta} & :=& \underset{\theta \in \Theta}{\rm argmin} ~ L(\theta) \tag{7} \\ \\ D_{KL}( \theta_{0} \mid\mid \hat{\theta} ) & := & \int f(x|\theta_{0}) \log \frac{f(x|\theta_{0})}{f(x|\hat{\theta})} ~dx \\ & = & \mathbb{E}_{X} \left[ \log \frac{ f(x|\theta_{0}) }{ f(x|\hat{\theta}) } \right] \tag{8.1} \\ & = & \mathbb{E}_{X} \bigl[ \ell( \theta_{0}|x ) \bigr] - \mathbb{E}_{X} \bigl[ \ell( \hat{\theta}|x ) \bigr] \tag{8.2} \end{eqnarray}

2.3 KL-divergenceとFisher情報量

　確率モデル $f(\cdot | \theta)$ を仮定したとき、パラメータ $\theta$ のFisher情報量 $I(\theta)$ は以下のように定義されます。（ただし、 $\ell( \cdot | x )$ を対数尤度関数とする。）

\begin{eqnarray} I(\theta) & := & \mathbb{E}_{X} \left[ { \left\{ \frac{d}{dx}\ell(\theta | x) \right\} }^{2} \right] \tag{9.1} \\ & = & \mathbb{E}_{X} \left[ { \left\{ \frac{d}{dx}\log f(x|\theta) \right\} }^{2} \right] \tag{9.2} \end{eqnarray}

　　確率モデル $f(\cdot | \theta)$ を仮定したとき、KL-divergenceとFisher情報量 $I(\theta)$ の間には、次のような関係*5が成り立ちます。

\begin{eqnarray} \lim_{h \to 0} \frac{1}{h^{2}} D_{KL} \bigl( f(x|\theta) \mid\mid f(x|\theta+h) \bigr) = \frac{1}{2} I(\theta) \tag{10} \end{eqnarray}

　この式は、 $\theta$ の近傍でのKL-divergence： $D_{KL} \bigl( f(x|\theta) \mid\mid f(x|\theta+h) \bigr)$ は、 $\theta$ における確率モデル $f(\cdot|\theta)$ のFisher情報量 $I(\theta)$ に比例することを示しています。つまり、Fisher情報量 $I(\theta)$ は確率モデル $f(\cdot | \theta)$ が $\theta$ の近傍でもつ局所的な情報量を表していることがわかります。

つまり，パラメトリックな確率モデル $p(x \vert \theta)$ に対するKL-divergenceの微小変化はFisher情報行列 $I(\theta)$ を計量行列として計算できます．

\begin{eqnarray} D_{KL}(p_{\theta}, p_{\theta + \Delta \theta}) &= \int p(x|\theta) \log \frac{p(x|\theta)}{p(x|\theta+\Delta \theta)} dx \ &\approx \frac{1}{2} {\Delta \theta}^T I(\theta) \Delta \theta \end{eqnarray}

f:id:yumaloop:20190107190651p:plain:w600

3. 参考書籍

新装改訂版現代数理統計学

作者:彰通, 竹村
発売日: 2020/11/10
メディア: 単行本

現代数理統計学の基礎 (共立講座数学の魅力)

作者:達也, 久保川
発売日: 2017/04/07
メディア: 単行本

数理統計学 (数学シリーズ)

作者:宣生, 稲垣
発売日: 2003/02/25
メディア: 単行本

Elements of Information Theory (Wiley Series in Telecommunications and Signal Processing)

作者:Cover, Thomas M.,Thomas, Joy A.
発売日: 2006/06/30
メディア: ハードカバー

情報量規準 (シリーズ・予測と発見の科学)

作者:貞則, 小西,源四郎, 北川
発売日: 2004/09/01
メディア: 単行本

*1:一般化するとf-divergenceというクラスが定義できます。汎関数微分などの関数解析の知識が必要。。

*2:pythonのライブラリ関数scipy.stats.entropyを使いました。

*3:熱力学的なエントロピーはボルツマン発祥ですが、シャノンの情報量に関する歴史的背景については以下で言及されています。ハートレー→ナイキスト→シャノンという流れがあるみたい。http://www.ieice.org/jpn/books/kaishikiji/200112/200112-9.html

*4:情報量規準についての文献→ https://www.ism.ac.jp/editsec/toukei/pdf/47-2-375.pdf

*5:この式は、KL-divergenceをテイラー展開して二次形式で近似すると導出されます。$$ \ell(\theta + h) - \ell(\theta) = {\ell}^{'}(\theta)h + \frac{1}{2} {\ell}^{''}(\theta) h^{2} + O(h^{3}) $$

2018-12-25

最尤推定・MAP推定・ベイズ推定を比較する

統計学

1. 推定のモチベーション
2. 最尤推定（ML）
- 2.1 パラメータの対数尤度
- 2.2 パラメータの推定（最尤推定量）
3. 事後確率最大化推定（MAP）
- 3.1 パラメータの事後確率
- 3.2 パラメータの推定（MAP推定量）
4. ベイズ推定（Bayse）
- 4.1 パラメータの平均対数尤度
- 4.2 確率分布の推定（予測分布）
5. 参考図書

1. 推定のモチベーション

　確率変数 $X$ の真の分布 $q(X)$ に対して、

$$ \begin{cases} 確率モデル　：　p(X|w) \\ 事前分布　　：　\varphi(w) \end{cases} $$

を仮定したとき、

$$ データサンプル　　：　X^{n} := \left\{ X_1, \dots ,X_n \right\} $$

を用いて、 $q(X)$ を推定したい。

2. 最尤推定（ML）

2.1 パラメータ $w$ の対数尤度 $\ell(w|\cdot)$

ベイズの定理： $$ \begin{eqnarray} p(w | X^{n}) = \frac{ p(X^{n}|w) \varphi(w) }{ \int_{W} p(X^{n}|w) \varphi(w) dw} \end{eqnarray} $$

において、以下で定義される対数尤度 $\ell(w|\cdot)$ は、手元にあるサンプルデータ $X^{n}$ に対する、確率モデル $p(x|w)$ の尤もらしさを示します。対数尤度 $\ell(w|\cdot)$ の示す値は、サンプルデータ $X^{n}$ がもつ"ばらつき"が加味されておらず、手元にあるデータの値に大きく依存するということに注意してください。

$$ \begin{eqnarray} \ell(w | \cdot) & := & \log p(\cdot | w) \\ \\ \ell(w | X^{n}) & := & \log p(X^{n} | w) \\ & = & \log \prod_{i=1}^{n} p(X_i | w) \end{eqnarray} $$

2.2 パラメータの推定（最尤推定量）

　最尤推定では、対数尤度 $\ell$ が最大となるようなパラメータ $\hat{w}_{ML}$ を選び、これを確率モデル $p(X|w)$ におけるパラメータ $w$ の推定量とします。推定量 $\hat{w}_{ML}$ を最尤推定量と呼びます。

$$ \begin{eqnarray} \hat{w}_{ML} := \underset{w \in W}{\rm argmax} ~ \ell(w|X^{n}) \end{eqnarray} $$

以上から、最尤推定では、確率変数 $X$ の真の分布 $q(x)$ は、サンプルデータ $X^{n}$ から求めた最尤推定量 $\hat{w}_{ML}$ と確率モデル $p(X|w)$ によって次式のように推定されます。

$$ \begin{eqnarray} q(X) \approx p( X | \hat{w}_{ML} ) \end{eqnarray} $$

3. 事後確率最大化推定（MAP）

3.1 パラメータ $w$ の事後確率 $p(w|\cdot)$

ベイズの定理：

$$ \begin{eqnarray} p(w | X^{n}) = \frac{ p(X^{n}|w) \varphi(w) }{ \int_{W} p(X^{n}|w) \varphi(w) dw} \end{eqnarray} $$

において、以下で定義される事後確率 $p(w|\cdot)$ は、手元にあるサンプルデータ $X^{n}$ に対する、確率モデル $p(X|w)$ ・事前分布 $\varphi(w)$ の尤もらしさを示します。事後確率 $p(w|\cdot)$ の示す値は、サンプルデータ $X^{n}$ がもつ"ばらつき"が加味されておらず、その値に大きく依存するということに注意してください。

$$ \begin{eqnarray} p(w | \cdot) & := & \frac{1}{Z} \log p(\cdot | w)\varphi(w) \\ \\ p(w | X^{n}) & := & \frac{1}{Z_n} \log p(X^{n} | w)\varphi(w) \\ & = & \frac{1}{Z_n} \log \prod_{i=1}^{n} p(X_i | w)\varphi(w) \end{eqnarray} $$

ただし、 $Z,~Z_n$ は、それぞれベイズの定理における右辺の分子に対して $w$ で積分した値（正規化定数）を表す。

3.2 パラメータの推定（MAP推定量）

　MAP推定では、事後確率が最大となるようなパラメータ $\hat{w}_{MAP}$ を選び、これを確率モデル $p(X|w)$ におけるパラメータ $w$ の推定量とします。推定量 $\hat{w}_{MAP}$ をMAP推定量と呼びます。

$$ \begin{eqnarray} \hat{w}_{MAP} := \underset{w \in W}{\rm argmax} ~ p(w|X^{n}) \end{eqnarray} $$

　以上から、MAP推定では、確率変数 $X$ の真の分布 $q(x)$ は、サンプルデータ $X^{n}$ から求めたMAP推定量 $\hat{w}_{ML}$ と確率モデル $p(X|w)$ によって次式のように推定されます。
$$ \begin{eqnarray} q(X) \approx p( X | \hat{w}_{MAP} ) \end{eqnarray} $$

4. ベイズ推定（Bayse）

4.1 パラメータ $w$ の平均対数尤度 $K$

ベイズの定理： $$ \begin{eqnarray} p(w | X^{n}) = \frac{ p(X^{n}|w) \varphi(w) }{ \int_{W} p(X^{n}|w) \varphi(w) dw} \end{eqnarray} $$

において、以下で定義される平均対数尤度 $K$ は、 $X$ に対して期待値をとった（ $X$ のばらつきによる偶然誤差を排除した）上での、確率モデル $p(X|w)$ の尤もらしさを示します。

$$ \begin{eqnarray} K(w) & := & \mathbb{E}_{X} \left[ \log p(X | w) \right] \tag{期待対数尤度} \\ \\ K_n(w) & := & \frac{1}{n} \sum_{i=1}^{n} \log p(X_i | w) \tag{経験対数尤度} \end{eqnarray} $$

パラメータの推定量を求めるため、平均対数尤度を符号反転したものを、損失関数 $L$ として定義します。

$$ \begin{eqnarray} L(w) & := & - K(w) = \mathbb{E}_{X} \left[ - \log p(X | w) \right] \tag{期待損失関数} \\ \\ L_n(w) & := & - K_{n}(w) = \frac{1}{n} \sum_{i=1}^{n} \left\{ - \log p(X_i | w) \tag{経験損失関数} \right\} \end{eqnarray} $$

4.2 確率分布 $q(x)$ の推定（予測分布）

　ベイズ推定では、損失関数が最小となるようなパラメータ $\hat{w}_{0}$ を選び、これを確率モデル $p(X|w)$ と事前分布 $\varphi(w)$ によって求められる事後分布 $p(w|X^{n})$ の平均値の推定量とみなします。ここで、最尤推定やMAP推定では、「確率モデル $p(X|w)$ に対して与えるパラメータ $w$ はある1つの値に決まる」と仮定していますが、ベイズ推定では「パラメータ $w$ は事前分布 $\varphi(w)$ によって確率的に変動する値（=確率変数）である」と仮定していることに注意します。推定量 $\hat{w}_{0}$ をベイズ推定量 *1と呼びます。

$$ \begin{eqnarray} \hat{w}_{0} := \underset{w \in W}{\rm argmax} ~ L_{n}(w) \end{eqnarray} $$

　推定量 $\hat{ w }_{0}$ は「事後分布において $w$ がとり得ると期待される（尤もらしい）値」と考えられます。しかし、ベイズ推定において、確率変数 $X$ の真の分布 $q(x)$ を、 $p( X | \hat{w}_{MAP} )$ によって推定することはしません。その代わり、推定量 $\hat{ w }_{0}$ を平均にもつ事後分布 $p(w|X^{n})$ から求められる予測分布 $p^{*}(X)$ を用いて真の分布 $q(x)$ を推定します。

$$ \begin{eqnarray} q(X) \approx p^{*}(X) & := & \int_{W} p(X|w)p(w|X^{n}) dw \end{eqnarray} $$

　上図は、確率モデル・事前分布・事後分布がすべて正規分布に従う場合の各分布を表す。Prior（事前分布）の「ばらつき」に対して、Likelihood（尤度, 最尤推定で用いる）と Posterior（事後分布, MAP推定で用いる）は、より小さい「ばらつき」をもつ。

5. 参考図書

StanとRでベイズ統計モデリング (Wonderful R)

作者: 松浦健太郎,石田基広
出版社/メーカー: 共立出版
発売日: 2016/10/25
メディア: 単行本
この商品を含むブログ (10件) を見る

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (29件) を見る

ベイズ統計の理論と方法

作者: 渡辺澄夫
出版社/メーカー: コロナ社
発売日: 2012/03/01
メディア: 単行本
購入: 1人クリック: 4回
この商品を含むブログ (8件) を見る

*1:他の言い方もあるようです。

2018-12-22

確率変数に関する4つの「収束概念」

統計学

1. 確率収束 : convergence in probability
2. p次平均収束 : convergence in mean
3. 分布収束（法則収束）: convergence in distribution (law)
4. 概収束 : almost sure convergence
参考文献

　確率変数の収束（Convergence of random variables）についてまとめてみました。数理統計学において極限現象の解析は特に重要で、漸近理論、極値理論、大偏差理論などで収束概念は欠かせません。
　解析学において、数列の極限から関数の極限を定義したように、統計学（確率論）においても確率変数列の極限から確率変数の極限を定義します。

1. 確率収束 : convergence in probability

確率変数列 $U_n := \left\{ U_1, U_2, \dots \right\}$ に対して、任意の(どんなに小さな) ${\varepsilon}>0$ について、 $$ \begin{eqnarray} \lim_{n \to \infty} P(~|U_n-U|{\geq}{\varepsilon}~) = 0 \end{eqnarray} $$ が成り立つとき、「確率変数列 $\left\{ U_n \right\}$ は、確率変数 $U$ に確率収束する」といい、これを $$ \begin{eqnarray} U_n \underset{~~~~~p}{\to} U \end{eqnarray} $$ と表記する。

2. p次平均収束 : convergence in mean

確率変数列 $U_n := \left\{ U_1, U_2, \dots \right\}$ に対して、 $$ \begin{eqnarray} \lim_{n \to \infty} E \left[~{|U_n-U|}^{~p} \right] = 0 \end{eqnarray} $$ が成り立つとき，「確率変数列 $\left\{ U_n \right\}$ は、確率変数 $U$ にp次平均収束する」という。　
　

3. 分布収束（法則収束）: convergence in distribution (law)

確率変数列 $U_n := \left\{ U_1, U_2, \dots \right\}$ に対して、確率変数 $U$ の分布関数 $F_U(x)$ の任意の連続点 $x$ で、 $$ \begin{eqnarray} \lim_{n \to \infty} P( U_n \leq x ) & = & P( U \leq x) \\ & = & F_U(x) \end{eqnarray} $$ が成り立つとき、「確率変数列 $\left\{ U_n \right\}$ は、確率変数 $U$ に分布収束（法則収束）する」といい、これを、 $$ \begin{eqnarray} U_n \underset{~~~~~d}{\to} U \end{eqnarray} $$ と表記する。　
　

4. 概収束 : almost sure convergence

確率変数列 $U_n := \left\{ U_1, U_2, \dots \right\}$ に対して、 $$ \begin{eqnarray} P\left( \omega~{\large |}~{\lim_{n \to \infty} \left| U_n(\omega)-U( \omega ) \right| =0} \right) = 1 \end{eqnarray} $$

が成り立つとき，「確率変数列 $\left\{ U_n \right\}$ は、確率変数 $U$ に概収束する」といい，これを、 $$ \begin{eqnarray} U_n \to U ~~~ a.s. \end{eqnarray} $$ と表記する。　　　　　

参考文献

現代数理統計学の基礎 (共立講座数学の魅力)

作者: 久保川達也
出版社/メーカー: 共立出版
発売日: 2017/04/07
メディア: 単行本
この商品を含むブログ (1件) を見る

現代数理統計学 (創文社現代経済学選書)

作者: 竹村彰通
出版社/メーカー: 創文社
発売日: 1991/12/01
メディア: 単行本
購入: 2人クリック: 26回
この商品を含むブログ (24件) を見る

入門確率過程

作者: 松原望
出版社/メーカー: 東京図書
発売日: 2003/11
メディア: 単行本
購入: 11人クリック: 197回
この商品を含むブログ (25件) を見る

2018-12-21

ベイズ推定の流れをまとめる（尤もらしさと不確実性）

統計学

ベイズ推定とは？
1. 事後分布
2. 予測分布
3. 損失関数
4. 誤差関数
5. 汎化損失・経験損失

ベイズ推定とは？

ベイズ推定の目的と主眼をまとめると、次のようになります。

ベイズ推定の目的
確率変数 $X$ の真の分布 $q(x)$ を、予測分布 $p^{*}(x)$ で推定する。
また、推定プロセスで生じる誤差や、推定結果に対する信頼性を評価したい。

このポストでは、『ベイズ統計の理論と方法』（渡辺澄夫, 2012）の内容を基に、ベイズ推定の基本の流れを整理するとともに、ベイズ推定で行われる「尤度最大化」と「情報量最小化」は同じであることを軸に、各統計量についての解説・解釈を述べたいと思います。

ベイズ統計の理論と方法

作者: 渡辺澄夫
出版社/メーカー: コロナ社
発売日: 2012/03/01
メディア: 単行本
購入: 1人クリック: 4回
この商品を含むブログ (8件) を見る

これから言及する「ベイズ推定の流れ」を、簡単にまとめました。

「ベイズ推定の流れ」

真の分布　　　:　 $q(x)$

に対して、

確率モデル　　：　 $p(x|w)$

事前分布　　　：　 $\varphi(w)$

を仮定したとき、

データ（標本） $X^{n} := \left\{ X_1, \dots ,X_n \right\}$

を用いて、以下の分布・統計量を求める。

事後分布　　　：　 $p(w|X_n)$

予測分布　　　：　 $p^{*}(x)$

期待損失関数　：　 $L(w)$

経験損失関数　：　 $L_n (w)$

期待誤差関数　：　 $K(w)$

経験誤差関数　：　 $K_n(w)$

汎化損失　　　：　 $G_{n}$

経験損失　　　：　 $T_{n}$

1. 事後分布 $p(w|X_n)$

　まず、確率モデル $p(x|w)$ と事前分布 $\varphi(x)$ から、パラメータ $w$ の事後分布 $p(w|X_n)$ を導出します。ベイズの定理より、

$$ \begin{eqnarray} p(w|X_n) & = & \frac{ p(X_{n}|w)\varphi(w) }{ p(X_{n}) } \\ \\ p(w|X_n) & = & \frac{ p(X_{n}|w)\varphi(w) }{ \int_{w} p(X_{n}|w)\varphi(w) dw } \end{eqnarray} $$

となります。さらに、右辺の分母を正規化定数 $Z$ とおけば、事後分布は以下のように定義することができます。

$$ \begin{eqnarray} p(w|X_n):=\frac{1}{Z} p(X_n|w)\varphi(w) \end{eqnarray} $$ ただし、 $$ \begin{eqnarray} Z:=\int_{w} p(X_{n}|w)\varphi(w) dw \end{eqnarray} $$

2. 予測分布 $p^{*}(x)$

　事前分布 $p(w)$ を事後分布 $p(X_n)$ におきかえて、再び確率モデル $p(X_n|w)$ のパラメータ $w$ に対する期待値をとります。

$$ \begin{eqnarray} p^{*}(x)=p(x|X_n):=\int_{w} p(x|w)p(w|X_n) dw \end{eqnarray} $$

ここまでの流れをまとめると、次のようになります。

確率モデル $p(x|w)$ と事前分布 $\varphi(w)$ を設定する。
確率モデル $p(x|w)$ と事前分布 $\varphi(w)$ によって、事後分布 $p(w|X_{n})$ を導出する。
確率モデル $p(x|w)$ と事後分布 $p(w|X_{n})$ によって、予測分布 $p^{*}(x)$ を導出する。
予測分布 $p^{*}(x)$ によって、真の分布 $q(x)$ を推定する。

次に、確率モデルの推定誤差を定量化する統計量（標本に対する関数） $L(~), K(~)$ や、推定誤差を小さくするパラメータ $w$ の推定量 $w_0$ を求めてみましょう。

3. 損失関数 $L(~)$

　「パラメータ $w$ の対数尤度の $X$ に対する期待値」として、損失関数 $L(~)$ を定義します。損失関数 $L(~)$ は、確率モデル $p(x|w)$ とパラメータ $w$ によって、決定されます。 $L(w)$ を期待損失関数、 $L_n(w)$ を経験損失関数と呼びます。実際には、真の分布 $q(x)$ がわからず期待損失関数 $L(w)$ を求められないので、データ $X^{n}$ を用いて、経験損失関数 $L_n(w)$ を導出することになります。

$$ \begin{eqnarray} L(w) & := & E_{X} \left[ -\log p(X|w) \right] \\ \\ L_n(w) & := & \frac{1}{n} \sum_{i=1}^{n} -\log p(X_i|w) \end{eqnarray} $$

よって、データ $X^{n}$ から導出される経験損失関数 $L_n(w)$ に対して、パラメータ $w$ の最適化問題を考えることができます。確率モデル $p(x|w)$ における、パラメータ $w \in W$ の推定量 $w_0$ を以下で定めます。

$$ \begin{eqnarray} w_0:=\underset{w \in W}{\rm argmin} ~ L_n(w) \end{eqnarray} $$

推定量 $w_0$ について：
　 $-\log p(X_i|w)$ は、確率モデル $p(X_i|w)$ によって推定された「データ $X_i$ に対する情報量（曖昧さ、不確実性）」を示しています。すなわち、ベイズ推定における推定量 $w_0$ は、「データ $X^{n}$ に対する不確実性を、平均的に最小化するようなパラメータ $w$ 」であると意味づけられます。
損失関数 $L(w)$ について：
　 $L(w)$ は、「あるパラメータ $w$ を与えたとき、確率モデル $p(X^{n}|w)$ がデータ $X^{n}$ に対してもつ平均的な不確実性」を表していると解釈できます。

4. 誤差関数 $K(~)$

「任意のパラメータ $w \in W$ に対する推定量 $w_0$ との誤差」として誤差関数 $K(~)$ を定義します。まず、パラメータ $w$ とその推定量 $w_0$ との誤差を、次の対数尤度比で定義します。

$$ \begin{eqnarray} f(x,w):= \log \frac{ p(x|w_0) }{ p(x|w) } \end{eqnarray} $$

　対数尤度比 $f(x, w)$ を用いて、誤差関数 $K(~)$ を定義します。 $K(w)$ を期待誤差関数、 $K_n(w)$ を経験誤差関数と呼びます。ここで、実際には真の分布 $q(x)$ がわからず期待誤差関数 $K(w)$ を求められないから、データ $X^{n}$ を用いて、経験誤差関数 $K_n(w)$ を導出することになります。

$$ \begin{eqnarray} K(w) & := & E_{X} \left[ f(X,w) \right]=E_{X} \left[ \log \frac{ p(X|w_0) }{ p(X|w) } \right] \\ \\ K_{n}(w) & := & \frac{1}{n} \sum_{i=1}^{n} f(X_i, w)=\frac{1}{n} \sum_{i=1}^{n} \log \frac{p(X_i|w_0)}{p(X_i|w)} \end{eqnarray} $$

対数尤度比 $f(x, w)$ について：
対数尤度比を分解すると、 $$ \begin{eqnarray} f(x, w) := \log \frac{ p(x|w_0) }{ p(x|w) } = \left[ -\log p(x|w) \right] - \left[ -\log p(x|w_0) \right] \geq 0 \end{eqnarray} $$ と表されるため、これは、あるパラメータ $w$ を与えたときの確率モデル $p(x|w)$ がもつ「尤もらしさ」を、「尤もらしさの上限」で正規化した値を示します。別の見方をすると、対数尤度比： $f(x, w)$ は、確率モデル $p(X^{n}|w)$ がもつ「データ $X^{n}$ に対する不確実性」を、「最小の不確実性」で正規化した値とも考えられます。
誤差関数 $K(w)$ について：
$K(w)$ は、「あるパラメータ $w$ を与えたとき、確率モデル $p(X^{n}|w)$ がデータ $X^{n}$ に対してもつ平均的な（正規化された）尤もらしさ」を表していると解釈できます。

5. 汎化損失・経験損失

　「予測分布 $p^{*}(x)$ と真の分布 $q(x)$ との誤差」として、「推定方法」がもつ損失を定義します。 $G_n$ を汎化損失、 $T_n$ を経験損失と呼びます。

$$ \begin{eqnarray} G_{n} & := & E_{X} \left[ -\log p^{*}(X) \right] \\ \\ T_{n} & := & \frac{1}{n} \sum_{i=1}^{n} - \log p^{*}(X_i)　 \end{eqnarray} $$

汎化損失 $G_n$ について：
$-\log p^{*}(X_i)$ は、予測分布 $p^{*}(X)$ によって推定された「データ $X_i$ に対する情報量（曖昧さ、不確実性）」を示します。すなわち、汎化損失 $G_n$ は、「パラメータ $w$ とデータ $X_i$ を与えたとき、確率モデル $p(X_n|w)$ から推定された予測分布 $p^{*}(X)$ が、データ $X^{n}$ に対してもつ平均的な不確実性」を表していると解釈できます。
汎化損失 $G_n$ と KL-距離 $D_{KL}(q \parallel p^{*})$ の比例関係
$$ \begin{eqnarray} G_n &=& E_{X} [ -\log p^{*}(X) ] \\ \\ &=& \int_{X} q(x) \cdot -\log p^{*}(X) dx \\ \\ &=& - \int_{X} q(x) \log q(x) dx + \int_{X} q(x) \log \frac{q(x)}{p^{*}(x)} dx \\ \\ &=& H(X)+D_{KL}(q \parallel p^{*}) \end{eqnarray} $$ ここで、確率モデル $p(x|w)$ と上式の関係を考えると、 $H(X)$ は、確率モデル $p(x|w)$ に依らないから無視できます。つまり、汎化損失 $G_n$ は $X$ の真の分布 $q(x)$ と予測分布 $p^{*}(x)$ の間のカルバック・ライブラー距離（Kullback–Leibler divergence） $D_{KL}(q \parallel p^{*})$ は、確率モデル $p(x|w)$ に対して比例関係にあります。

$$ \begin{align} G_n \propto D_{KL}(q \parallel p^{*}) \end{align} $$

KL-距離 $D_{KL}(q \parallel p^{*})$ について：
真の分布 $q(x)$ と予測分布 $p^{*}(x)$ の間のカルバック・ライブラー距離（Kullback–Leibler divergence）の定義式は、 $$ \begin{eqnarray} D_{KL}(q \parallel p^{*}) := \int_{X} q(x) \log \frac{q(x)}{p^{*}(x)} dx \end{eqnarray} $$ ですから、これは「対数尤度比 $\log \frac{q(x)}{p^{*}(x)}$ の確率変数 $X$ に対する期待値」として解釈できます。例として、下図右は、2つの正規分布 $p(x)$ と $q(x)$ について、それぞれの分布に従う確率変数 $X$ の各値を横軸に、それによって計算される対数尤度比を縦軸にとったグラフを示しています。よって、 $D_{KL}(q \parallel p)$ は、その期待値ですから、下図右で青く塗りつぶされた部分の「確率重みつき面積（期待値）」となります。

2018-12-01

統計的仮説検定をていねいに解説する

統計学 CheatSheet 備忘

1. 推測統計学の概念・用語
- 1.1.　推測統計学とそのモチベーション
- 1.2.　概念・用語の定義
2. 統計的仮説検定
3. 正規母集団についての仮説検定

1. 推測統計学の概念・用語

仮説検定は、推測統計学の１分野です。そこで、仮説検定の説明に入る前に、まず推測統計学に関する概念や用語のかんたんな定義を述べます。

1.1.　推測統計学とそのモチベーション

　推測統計学は、記述統計学を発展させた統計学の体系で、現代的な統計学の基本です。推測統計学の土台には、大数の法則と中心極限定理があり、これらの定理によって分布の収束や近似といった極限的なふるまいの解析学的な扱いが可能になります。

　記述統計学のモチベーションは、「観察や調査によって集められた大量のデータそのものを目的に応じて整理・要約する」というものですが、推測統計学では「データの発生元である母集団に確率モデルを想定し、データを標本として捉えることで、データ自身ではなくデータを発生させている根源的な事象に対して推論を行う」という、より貪欲なモチベーションを置いています。

1.2.　概念・用語の定義

　手元にあるデータに対して、「母集団と標本」という "枠組み"（frame work）を与えることで、手元のデータと確率モデルを結びつけ、以降の議論を進めやすくします。

母集団と標本

母集団（population）
ーある確率分布、確率モデルに従うであろう事象（event）の全集合*1

標本（sample）
ー母集団の部分集合。母集団から抽出された確率変数（random variable）の組

データ（data）
ー標本の実現値（realized value）*2

　さらに、母集団や標本をそのまま捉えようとするのではなく、母集団の性質を示す「母数」や標本の性質を示す「統計量」を考えることで、理論の厳密さを上げます。（機械学習の文脈では、統計量のことを特徴量と呼ぶ。）

母数（parameter）
ー母集団が持っている特性を示す数値
　確率変数に依存しない定数 $\theta$ として定義される。*3

統計量（statistics）
ー標本が持っている特性を示す数値
　標本についての関数 $T(X_1, ... X_n)$ として定義される。統計量は確率変数。

次に、統計的推測をまとめ、推測統計学における「仮説検定」の位置付けをみてみましょう。統計的推測の手法には、目的に応じて「点推定」「区間推定」「仮説検定」「予測」などがあります。

統計的推測のまとめ：

統計的推測（statistical inference）
ー標本 $X_1, ... X_n$ に基づいて母集団（の確率分布や、確率モデル）に関する推論を行うこと。

点推定（estimation）
ー母数 $\theta$ の値を、標本 $X_1, ... X_n$ に基づいて言い当てること

区間推定（interval estimation）
ー母数 $\theta$ が入っているような信頼性の高い区間を、標本 $X_1, ... X_n$ に基づいて言い当てること

仮説検定（hypothesis testing）
ー母集団に関する仮説（：母数 $\theta$ や確率分布 $f(\cdot~|~\theta)$ に対する記述）が正しいか否かを、標本 $X_1, ... X_n$ に基づいて判定すること

予測（prediction）
ー母集団に関する確率変数 $X$ の新しい実現値 $x$ を、標本 $X_1, ... X_n$ に基づいて推測すること。

　このように、仮説検定とは、母集団の特性（値）を定量的に求めようとする「推定」や「予測」とは異なり、母集団の特性に対する仮説が正しいことを「支持する or 支持しない」の２者択一で判定するための方法論です。しかるに、有意水準 $\alpha$ や $P$ 値といった指標は、この二者択一の判断に対して妥当性（信頼性）を保証するための数値でしかありません。

　どんな検定方式を使おうが、「仮説検定」の枠組みに入る統計手法である限り、データから得られる結果は、仮説を「支持する or 支持しない」のいずれかであることに注意してください。仮説検定では、判断結果（棄却or受容）に対する妥当性（信頼性）を見積もることができますが、データに対する仮説そのものの妥当性（信頼性）を見積もることはできません。

2. 統計的仮説検定

2.1.　母数と標本

　仮説検定で重要な概念となる「母数、母数空間」と「標本、標本空間」について、改めておさらいします。

母数空間 $\Theta$ と標本空間 $\mathcal{X}$ ：

母数（parameter） $\theta$
ー母集団が持っている特性を示す数値.

標本（sample） $X$
ー母集団の部分集合. 母集団から抽出された確率変数（random variable）の組.

母数空間（parameter space） $\Theta$
ー母数 $\theta$ のとり得る値の集合.

標本空間（sample space） $\mathcal{X}$
ー標本 $X$ のとり得る値の集合.

2.2.　仮説と検定

　ここまでで準備が整ったので、いよいよ統計的仮説検定について言及していきます。まずは仮説検定を「仮説」と「検定」に分けて、定義してみます。

仮説検定：

仮説（hypothesis）
ー母集団に関する何らかの記述。母数 $\theta$ が、母数空間 $\Theta$ を分割する2つの排反する部分集合 $\Theta_0$ , $\Theta_1$ のどちらに属するかを記述する。

帰無仮説（null hypothesis）　　　 $H_0$ : $\theta \in \Theta_0$ 「母数 $\theta$ が $\Theta_0$ に属する」

対立仮説（alternative hypothesis） $H_1$ : $\theta \in \Theta_1$ 「母数 $\theta$ が $\Theta_1$ に属する」

検定（testing）
ー帰無仮説 $H_0$ が「真か偽か」の判定を、標本 $X$ から求める検定統計量 $T(X)$ で行う。

受容（accept）：「帰無仮説 $H_0$ が真である」と検定する。

棄却（reject）：「帰無仮説 $H_0$ が偽である」と検定する。

ここで注意すべきは、以下の点です。

「仮説」は、母数空間 $\Theta$ に対する記述。
「検定」は、標本空間 $\mathcal{X}$ に基づいて行われる。

　さて、次に問題となるのは、「母集団に関する何らかの記述が与えられたとき、それをどのように「真 / 偽」と判定するか？」、言い換えると「ある仮説が与えられたとき、それをどのように検定するのか？」です。このような「検定を実施する方法」のことを、検定方式と言います。

　ある仮説に対する統計的仮説検定そのものの「性質」は、検定を行う方式（＝検定方式）によって決まります。１つの仮説（: 例えば「母集団の平均は1である」）に対して、「支持する or 支持しない」の判断を下す方法として、複数の検定方式が考えられます。

2.3.　有意水準 $\alpha$ , P値

母数空間 $\Theta$ において、母数 $\theta \in \Theta$ に関する仮説の検定を考える。

帰無仮説（null hypothesis）　　　 $H_0$ : $\theta \in \Theta_0$
対立仮説（alternative hypothesis） $H_1$ : $\theta \in \Theta_1$

A） 棄却域
帰無仮説 $H_0$ の棄却域 $R$ は、標本 $X \in \mathcal{X}$ に関して定められ、以下のように定義されます。

　　　 $R := {\large\{} {\it x} \in \mathcal{X} ~{\large |}~ \text{帰無仮説}~{H_0}~\text{が棄却されるための標本統計量}~T(X)~\text{の条件} {\large\}}$

ここで、標本 $X$ は「母数空間 $\Theta$ に基づく確率分布に従う母集団」から抽出されたものだから、その確率は $P_{\theta}$ に従うため、「標本 $X$ が帰無仮説 $H_0$ の棄却域 $R$ に属する確率」は、以下のように定義されます。

　　　標本 $X$ が帰無仮説 $H_0$ の棄却域 $R$ に属する確率＝ $P_{\theta} \left( X \in R \right)$

確率 [tex: P{\theta} \left( X \in R \right)] によって、「帰無仮説 $H_0$ を棄却するか否か」を決定します。このとき「帰無仮説 $H_0$ を棄却する」と判断する際の確率 [tex: P{\theta} \left( X \in R \right)] の最大値を「有意水準 $\alpha$ 」と呼びます。

有意水準 $\alpha$ の検定：

「標本 $X$ に基づいて、帰無仮説 $H_0$ : $\theta \in \Theta_0$ を棄却する」という判断が正しくない確率は、大きくても $\alpha$ 以下である。

$\underset{\theta \in \Theta_0}\sup P_{\theta} \left( X \in R \right) \leq \alpha$
　　

$H_0$ の棄却域 $R := {\large\{} {\it x} \in \mathcal{X} ~{\large |}~帰無仮説~{H_0}~が棄却されるための標本統計量~T(X)~の条件{\large\}}$

さらに標本 $X$ の棄却域 $R$ を書き下すと、

$\underset{\theta \in \Theta_0}\sup P_{\theta} \left( X \in {{\large\{} {\it x} \in \mathcal{X} ~{\large |}~ 帰無仮説~{H_0}~が棄却されるための標本統計量~T(X)~の条件 {\large\}}} \right) \leq \alpha$

となる。

f:id:yumaloop:20181202011259p:plain

3. 正規母集団についての仮説検定

3.1.　1標本の正規母集団

　 $X_1, \ldots, X_n ~~~ i.i.d. ~{\sim}~~ N( \mu,~{{\sigma}^2})$ とする。
　また、標本をまとめて $\mathbf{X}~=~\left( X_1, \ldots, X_n \right)$ とおく。

標本平均： $\overline{X}~=~{\large \frac{1}{n} } \sum_{i=1}^{n} X_i$
標本分散： ${V^{2}}~=~{\large \frac{1}{n-1} } \sum_{i=1}^{n} {\left( X_i - \overline{X} \right) }^2$

　以下のように、母平均 $\mu$ と母分散 ${\sigma}^2$ に関する仮説検定を考える。

1標本における, 母平均 $\mu$ と母分散 ${\sigma}^2$ についての仮説検定

(a) 母平均 $\mu$ の両側検定 $\begin{cases} \text{帰無仮説} ~ H_0: \mu = \mu_0 \\ \text{対立仮説} ~ H_1: \mu \neq \mu_0 \end{cases}$

(b) 母平均 $\mu$ の片側検定 $\begin{cases} \text{帰無仮説} ~ H_0: \mu = \mu_0 \\ \text{対立仮説} ~ H_1: \mu > \mu_0 \end{cases}$

(c) 母分散 ${\sigma}^2$ の両側検定 $\begin{cases} \text{帰無仮説} ~ H_0: {\sigma}^2 = {\sigma_0}^2 \\ \text{対立仮説} ~ H_1: {\sigma}^2 \neq {\sigma_0}^2 \end{cases}$

3.1.1.　母平均 $\mu$ の検定（母分散 ${\sigma}^{2}$ が既知）

　　　 $H_0: \mu = \mu_0$ の下で、中心極限定理より、

　　　 ${\large \frac{\overline{X}~-~{\mu_0}}{ \sqrt{\frac{{\sigma}^{2}}{n}} } }~ \overset{d}{\longrightarrow} ~N(0, 1)$

(a) 両側検定の帰無仮説 $H_0: \mu = \mu_0$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large |}~{\sqrt{n}}{\large \frac{ | \overline{X}~-~{\mu_0} | }{\sigma} } > {Z_{\frac{\alpha}{2}}} \right\}$

(b) 片側検定の帰無仮説 $H_0: \mu = \mu_0$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large |}~{\sqrt{n}}{\large \frac{ \overline{X}~-~{\mu_0} }{\sigma} } > {Z_{\alpha}} \right\}$

※ 母平均 $\mu$ の信頼係数 $1-\alpha$ の信頼区間

$C(\mathbf{X})=\left\{ \mu ~{\large |}~ {\overline{X}-{\large \frac{\sigma}{\sqrt{n}}}{Z_{\frac{\alpha}{2}}}} \leq ~ \mu ~ \leq {\overline{X}+{\large \frac{\sigma}{\sqrt{n}} } {Z_{\frac{\alpha}{2}}} } \right\}$

3.1.2.　母平均 $\mu$ の検定（母分散 ${\sigma}^{2}$ が未知）

　　　 $H_0: \mu = \mu_0$ の下で、

　　　 ${\large \frac{\overline{X}~-~{\mu_0}}{ \sqrt{\frac{{V}^{2}}{n}} } }~ {\longrightarrow} ~{t_{n-1}}$

(a) 両側検定の帰無仮説 $H_0: \mu = \mu_0$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large |}~{\sqrt{n}}{\large \frac{ | \overline{X}~-~{\mu_0} | }{V} } > {t_{{n-1},~{\frac{\alpha}{2}}}} \right\}$

(b) 片側検定の帰無仮説 $H_0: \mu = \mu_0$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large |}~{\sqrt{n}}{\large \frac{ \overline{X}~-~{\mu_0} }{V} } > {t_{{n-1},~{\alpha}}} \right\}$

※ 母平均 $\mu$ の信頼係数 $1-\alpha$ の信頼区間

$C(\mathbf{X})=\left\{ \mu ~{\large |}~ {\overline{X}-{\large \frac{V}{\sqrt{n}}}{t_{{n-1},~{\frac{\alpha}{2}}}}} ~{\leq} ~ {\mu} ~ {\leq}~ {\overline{X}+{\large \frac{V}{\sqrt{n}} } t_{n-1,~\frac{\alpha}{2}} } \right\}$

3.1.3.　母分散 ${\sigma}^2$ の検定

　　　 $H_0: \mu = \mu_0$ の下で、

　　　 ${\large \frac{(n-1){{V}^{2}}}{{{\sigma}_0}^2} } ~ {\sim} ~ \chi_{n-1}^{2}$

(c) 両側検定の帰無仮説 $H_0: {\sigma}^2 = {\sigma_0}^2$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large |}~{\large \frac{(n-1){{V}^2}}{{{\sigma}_0}^2} } \lt {\chi_{n-1,~1-{\frac{\alpha}{2}}}^{2}}, \ or \quad {\chi_{n-1,~1-{\frac{\alpha}{2}}}^{2}} \lt {\large \frac{(n-1){{V}^2}}{{{\sigma}_0}^2} } \right\}$

※ 母分散 ${\sigma}^2$ の信頼係数 $1-\alpha$ の信頼区間

$C(\mathbf{X})=\left\{ {{\sigma}^2} ~{\large |}~{\large \frac{(n-1){{V}^2}}{ \chi_{n-1,~{\frac{\alpha}{2}}}^2 }} \leq {{\sigma}^2} \leq {\large \frac{(n-1){{V}^2}}{ \chi_{n-1,~1-{\frac{\alpha}{2}}}^2 }} \right\}$

3.2　2標本の正規母集団（母分散が等しい）

　 $X_1, \ldots, X_m ~~~ i.i.d. ~{\sim}~~ N( \mu_1,~{{\sigma}^2})$ とする。
　 $Y_1, \ldots, Y_n ~~~ i.i.d. ~{\sim}~~ N( \mu_2,~{{\sigma}^2})$ とする。

　また、標本をまとめて、
　 $\mathbf{X}~=~\left( X_1, \ldots, X_m \right)$ ,
　 $\mathbf{Y}~=~\left( Y_1, \ldots, Y_n \right)$
　とおく。

$\mathbf{X}$ の標本平均： $\overline{X}~=~{\large \frac{1}{m} } \sum_{i=1}^{m} X_i$
$\mathbf{Y}$ の標本平均： $\overline{Y}~=~{\large \frac{1}{n} } \sum_{i=1}^{n} Y_i$

　以下のように、母平均 $\mu$ と母分散 ${\sigma}^2$ に関する仮説検定を考える。

2標本における, 母平均 $\mu_1,~\mu_2$ についての仮説検定

(a) 母平均 $\mu_1,~\mu_2$ の同等性についての両側検定 $\begin{cases} \text{帰無仮説} ~ H_0: \mu_1 = \mu_2 \\ \text{対立仮説} ~ H_1: \mu_1 \neq \mu_2 \end{cases}$

(b) 母平均 $\mu_1,~\mu_2$ の同等性についての片側検定 $\begin{cases} \text{帰無仮説} ~ H_0: \mu_1 = \mu_2 \\ \text{対立仮説} ~ H_1: \mu_1 > \mu_2 \end{cases}$

3.2.1　母平均 $\mu_1,~\mu_2$ の同等性の検定（母分散 ${{\sigma}^2}$ が既知）

　　　 $H_0: \mu_1 = \mu_2$ の下で、中心極限定理（CLT）より、

　　　 ${\large \frac{\overline{X}~-~\overline{Y}}{ \sqrt{\frac{m+n}{mn} {{\sigma}^2}} } }~ \overset{d}{\longrightarrow} ~N(0, 1)$

(a) 両側検定の帰無仮説 $H_0: \mu_1 = \mu_2$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large |}~{\large \sqrt{\frac{mn}{m+n}} \frac{ | \overline{X}~-~\overline{Y} | }{\sigma}} > {Z_{\frac{\alpha}{2}}} \right\}$

(b) 片側検定の帰無仮説 $H_0: \mu_1 = \mu_2$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large |}~{\large \sqrt{\frac{mn}{m+n}} \frac{ \overline{X}~-~\overline{Y} }{\sigma}} > {Z_{\alpha}} \right\}$

※ 母平均の差 $\mu_1 - \mu_2$ の信頼係数 $1-\alpha$ の信頼区間

$C(\mathbf{X})=\left\{ \mu_1 - \mu_2 ~{\large |}~ {(\overline{X} - \overline{Y})-{\sigma}{\large \sqrt{\frac{m+n}{mn}}}}{Z_{\frac{\alpha}{2}}} \lt {\mu_1 - \mu_2} \lt {(\overline{X} - \overline{Y})+{\sigma}{\large \sqrt{\frac{m+n}{mn}}}}{Z_{\frac{\alpha}{2}}} \right\}$

3.2.2　母平均 $\mu_1,~\mu_2$ の同等性の検定（母分散 ${{\sigma}^2}$ が未知）　"t検定"

　　　標本から（2標本に共通する）母分散を推定した値 ${{\hat{\sigma}}^2}$ を以下のように定義する。

　　　 ${{\hat{\sigma}}^2} := {\large \frac{1}{m+n-2}} \left\{ {\sum_{i=1}^{m} {\left( X_i - \overline{X} \right) }^2} + {\sum_{j=1}^{n} {\left( Y_j - \overline{Y} \right) }^2} \right\}$

　　　 $H_0: \mu_1 = \mu_2$ の下で、

　　　 ${\large \frac{\overline{X}~-~\overline{Y}}{ \sqrt{\frac{m+n}{mn} {{\hat{\sigma}}^2}} } }~ {\sim} ~{t_{m+n-2}}$

(a) 両側検定の帰無仮説 $H_0: \mu_1 = \mu_2$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large |}~{\large \sqrt{\frac{mn}{m+n}} \frac{ | \overline{X}~-~\overline{Y} | }{\sigma}} > {t_{{m+n-2},~ \frac{\alpha}{2}}} \right\}$

(b) 片側検定の帰無仮説 $H_0: \mu_1 = \mu_2$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large |}~{\large \sqrt{\frac{mn}{m+n}} \frac{ ( \overline{X}~-~\overline{Y} ) }{\sigma}} > {t_{{m+n-2},~ \alpha}} \right\}$

※ 母平均の差 $\mu_1 - \mu_2$ の信頼係数 $1-\alpha$ の信頼区間

$C(\mathbf{X})=\left\{ \mu_1 - \mu_2 ~{\large |}~ {(\overline{X} - \overline{Y})-{\sigma}{\large \sqrt{\frac{m+n}{mn}}}}{t_{m+n-2,~\frac{\alpha}{2}}} \lt {\mu_1 - \mu_2} \lt {(\overline{X} - \overline{Y})+{\sigma}{\large \sqrt{\frac{m+n}{mn}}}}{t_{m+n-2,~\frac{\alpha}{2}}} \right\}$

3.3　2標本の正規母集団（母分散が等しくない）

　 $X_1, \ldots, X_m ~~~ i.i.d. ~{\sim}~~ N( \mu_1,~{{\sigma_1}^2})$ とする。
　 $Y_1, \ldots, Y_n ~~~ i.i.d. ~{\sim}~~ N( \mu_2,~{{\sigma_2}^2})$ とする。

　また、標本をまとめて、
　 $\mathbf{X}~=~\left( X_1, \ldots, X_m \right)$ ,
　 $\mathbf{Y}~=~\left( Y_1, \ldots, Y_n \right)$
　とおく。

$\mathbf{X}$ の標本平均： $\overline{X}~=~{\large \frac{1}{m} } \sum_{i=1}^{m} X_i$
$\mathbf{Y}$ の標本平均： $\overline{Y}~=~{\large \frac{1}{n} } \sum_{i=1}^{n} Y_i$
$\mathbf{X}$ の標本分散： ${V_{1}^{2}}~=~{\large \frac{1}{m-1} } \sum_{i=1}^{m} {\left( X_i - \overline{X} \right) }^2$
$\mathbf{Y}$ の標本分散： ${V_{2}^{2}}~=~{\large \frac{1}{n-1} } \sum_{i=1}^{n} {\left( Y_i - \overline{Y} \right) }^2$

　以下のように、母平均 $\mu$ と母分散 ${\sigma}^2$ に関する仮説検定を考える。

2標本における, 母平均 $\mu$ と母分散 ${\sigma}^2$ についての仮説検定

(a) 母平均 $\mu_1,~\mu_2$ の同等性についての両側検定 $\begin{cases} \text{帰無仮説} ~ H_0: \mu_1 = \mu_2 \\ \text{対立仮説} ~ H_1: \mu_1 \neq \mu_2 \end{cases}$

(b) 母平均 $\mu_1,~\mu_2$ の同等性についての片側検定 $\begin{cases} \text{帰無仮説} ~ H_0: \mu_1 = \mu_2 \\ \text{対立仮説} ~ H_1: \mu_1 \gt \mu_2 \end{cases}$

(c) 母分散 ${\sigma_1}^2,~{\sigma_2}^2$ の同等性についての片側検定 $\begin{cases} \text{帰無仮説} ~ H_0: {\sigma_1}^2 = {\sigma_2}^2 \\ \text{対立仮説} ~ H_1: {\sigma_1}^2 \neq {\sigma_2}^2 \end{cases}$

3.3.1　母平均 $\mu_1,~\mu_2$ の同等性の検定（母分散 ${{\sigma_1}^2},~{{\sigma_2}^2}$ が既知）

　　　 $H_0: \mu_1 = \mu_2$ の下で、中心極限定理（CLT）より、

　　　 ${\large \frac{\overline{X}~-~\overline{Y}}{ \sqrt{\frac{{\sigma_1}^2}{m} + \frac{{\sigma_2}^2}{n}} } }~ \overset{d}{\longrightarrow} ~N(0, 1)$

(a) 両側検定の帰無仮説 $H_0: \mu_1 = \mu_2$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large |}~{\large \frac{ | \overline{X}~-~\overline{Y} | }{ \sqrt{\frac{{\sigma_1}^2}{m} + \frac{{\sigma_2}^2}{n}} }} > {Z_{ \frac{\alpha}{2}} } \right\}$

(b) 片側検定の帰無仮説 $H_0: \mu_1 = \mu_2$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large |}~{\large \frac{ ( \overline{X}~-~\overline{Y} ) }{ \sqrt{\frac{{\sigma_1}^2}{m} + \frac{{\sigma_2}^2}{n}} }} > {Z_{ \alpha }} \right\}$

※ 母平均の差 $\mu_1 - \mu_2$ の信頼係数 $1-\alpha$ の信頼区間

$C(\mathbf{X})=\left\{ \mu_1 - \mu_2 ~{\large |}~ {(\overline{X} - \overline{Y})-{\large \sqrt{\frac{{\sigma_1}^2}{m} {\small +} \frac{{\sigma_2}^2}{n}}}}{Z_{\frac{\alpha}{2}}} \lt {\mu_1 - \mu_2} \lt {(\overline{X} - \overline{Y})+{\large \sqrt{\frac{{\sigma_1}^2}{m} {\small +} \frac{{\sigma_2}^2}{n}}}}{Z_{\frac{\alpha}{2}}} \right\}$

3.3.2　母平均 $\mu_1,~\mu_2$ の同等性の検定（母分散 ${{\sigma_1}^2},~{{\sigma_2}^2}$ が未知）

　　　 $H_0: \mu_1 = \mu_2$ の下で、Welch's test（ウェルチの検定）を用いると、

　　　 ${\large \frac{\overline{X}~-~\overline{Y}}{ \sqrt{\frac{{V_1}^2}{m} + \frac{{V_2}^2}{n}} } }~ \underset{近似分布}{\longrightarrow} ~t_f$

　　　ただし、 $t$ 分布の自由度 $f$ は以下のように求める。

　　　 $f:={\large \frac{{( \frac{{V_1}^2}{m} + \frac{{V_2}^2}{n} )}^2}{ \frac{{V}_{2}^{4}}{m(m-1)}+\frac{{V}_{2}^{4}}{n(n-1)} } }$

(a) 両側検定の帰無仮説 $H_0: \mu_1 = \mu_2$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large \mid}~{\large \frac{ | \overline{X}~-~\overline{Y} | }{ \sqrt{\frac{{V_1}^2}{m} + \frac{{V_2}^2}{n}} }} > {t_{f,~ \frac{\alpha}{2}} } \right\}$

(b) 片側検定の帰無仮説 $H_0: \mu_1 = \mu_2$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large \mid}~{\large \frac{ ( \overline{X}~-~\overline{Y} ) }{ \sqrt{\frac{{V_1}^2}{m} + \frac{{V_2}^2}{n}} }} > {t_{f,~ \alpha} } \right\}$

※ 母平均の差 $\mu_1 - \mu_2$ の信頼係数 $1-\alpha$ の信頼区間

$C(\mathbf{X})=\left\{ \mu_1 - \mu_2 ~{\large \mid}~ {(\overline{X} - \overline{Y})-{\large \sqrt{\frac{{V_1}^2}{m} {\small +} \frac{{V_2}^2}{n}}}}{t_{f,~\frac{\alpha}{2}}} \lt {\mu_1 - \mu_2} \lt {(\overline{X} - \overline{Y})+{\large \sqrt{\frac{{V_1}^2}{m} {\small +} \frac{{V_2}^2}{n}}}}{t_{f,~ \frac{\alpha}{2}}} \right\}$

3.3.3　母分散 ${\sigma_1}^2,~{\sigma_2}^2$ の同等性の検定　"F検定"

　　　 $H_0: {\sigma_1}^2={\sigma_2}^2$ の下で、

　　　 ${\large \frac{ {V_1}^2 ~/~ {\sigma_1}^2 }{ {V_2}^2 ~/~ {\sigma_2}^2 } } = \frac{{V_1}^2}{{V_2}^2} ~ {\sim} ~F_{m-1,~n-1}$

(c) 両側検定の帰無仮説 $H_0: {\sigma_1}^2 = {\sigma_2}^2$ の有意水準 $\alpha$ の棄却域

$R=\left\{ \mathbf{X}\in\mathcal{X} ~{\large \mid}~ {\large \frac{{V_1}^2}{{V_2}^2}} \lt F_{m-1,~n-1,~1-\frac{\alpha}{2}} ,\ or \quad F_{m-1,~n-1,~1-\frac{\alpha}{2}} \lt {\large \frac{{V_1}^2}{{V_2}^2}} \right\}$

※ 母分散の比 ${\large \frac{{\sigma_1}^2}{{\sigma_2}^2}}$ の信頼係数 $1-\alpha$ の信頼区間

$C(\mathbf{X})=\left\{ {\large \frac{{\sigma_1}^2}{{\sigma_2}^2}} ~{\large \mid}~ {\frac{{V_1}^2}{{V_2}^2}} \frac{1}{F_{m-1,~n-1,~\frac{\alpha}{2}}} \leq {\large \frac{{\sigma_1}^2}{{\sigma_2}^2}} \leq {\frac{{V_1}^2}{{V_2}^2}} \frac{1}{F_{m-1,~n-1,~1-\frac{\alpha}{2}}} \right\}$

参考

現代数理統計学の基礎 (共立講座数学の魅力)

作者: 久保川達也
出版社/メーカー: 共立出版
発売日: 2017/04/07
メディア: 単行本
この商品を含むブログ (1件) を見る

統計学入門 (基礎統計学?)

作者: 東京大学教養学部統計学教室
出版社/メーカー: 東京大学出版会
発売日: 1991/07/09
メディア: 単行本
購入: 158人クリック: 3,604回
この商品を含むブログ (79件) を見る

*1:定義より、事象に対して、確率変数を定めることができる。厳密には、可測集合に含まれる事象に対して確率を定義し、確率に対して確率変数を定義する。

*2:「実現値」の概念は伝統的な数理統計学で扱うもので、ベイズ的な枠組みでは、あまり気にしない。例えば、以下のような批判もある。https://eprints.lib.hokudai.ac.jp/dspace/bitstream/2115/47712/1/ES_61(3)_001.pdf

*3:ベイズ統計学では、確率変数として扱う。

2018-08-27

Keras/Tensorflow : CIFAR-10のMobileNetV2-likeなアーキテクチャを作った

機械学習 Tensorflow Python Keras

f:id:yumaloop:20181201151138p:plain — MobilenetV2-残差ブロックの構造

0. 前提
1. 動作環境
2. プログラム
3. 学習結果

0. 前提

MobileNetV2は、2018/04にGoogleのMark Sandlerらによって発表された論文："MobileNetV2: Inverted Residuals and Linear Bottlenecks"にて導入された、ニューラルネットワークモデルです。その名の通り、モバイル機器向けの高速なフォワード処理に適しており、精度を保持した上でパラメータの削減に成功しています。これは、同じくGoogleから発表された論文："Mobilenets: Efficient Convolutional Neural Networks for Mobile Vision Applications"で提唱されたMobileNetアーキテクチャのアップグレード版という位置付けで、出力チャネル数を明に増加させる「Expand-layer」やResNet系で馴染み深い「skip-connection」の導入などが特徴的です。

MobileNetV2(原論文) : https://arxiv.org/pdf/1801.04381.pdf
MobileNetV2(解説Blog) : http://machinethink.net/blog/mobilenet-v2
Github - keras_applications :

github.com

1. 動作環境

OS: Ubuntu 16.04

Package             Version  
------------------- -------
python              3.5.0
tensorboard         1.9.0    
tensorflow          1.9.0    
h5py                2.8.0    
Keras               2.2.2    
Keras-Applications  1.0.4    
Keras-Preprocessing 1.0.2

2. プログラム

import_CIFAR-10.py

CIFAR-10データセットをnumpy配列（ndarray）形式で保存しておく実行ファイルです. CIFAR-10の画像は一枚あたり「32w(pixel) × 32h(pixel) × 3ch(RGB)」個のpixelからできています. 今回は, 画像データに対する前処理としてRGB-channelに応じて標準化（平均0, 標準偏差1）を行います.

データセット用のディレクトリと実行ファイルを用意しておくとバージョン管理が楽です.

import numpy as np
from keras.datasets import cifar10
from keras.utils import np_utils

nb_classes = 10
argmax_ch  = 255.0

if __name__=='__main__':
    # load CIFAR-10 data
    (X_train, y_train), (X_test, y_test) = cifar10.load_data()

    # set data type as 'float32'
    X_train = X_train.astype('float32') #argmax_ch
    X_test  = X_test.astype('float32')  #argmax_ch

    def ch_wise_normalization(X_type, ch):
        mean_ch = X_type[:, :, :, ch].mean()
        std_ch = X_type[:, :, :, ch].std()
        X_type[:, :, :, ch] = (X_type[:, :, :, ch] - mean_ch) / std_ch
        return X_type[:, :, :, ch]

    # normalize data for each R-G-B(0, 1, 2) channel 
    X_train[:, :, :, 0] = ch_wise_normalization(X_train, 0)
    X_train[:, :, :, 1] = ch_wise_normalization(X_train, 1)
    X_train[:, :, :, 2] = ch_wise_normalization(X_train, 2)

    X_test[:, :, :, 0]  = ch_wise_normalization(X_test, 0)
    X_test[:, :, :, 1]  = ch_wise_normalization(X_test, 1)
    X_test[:, :, :, 2]  = ch_wise_normalization(X_test, 2)

    # convert class label (0-9) to one-hot encoding format
    y_train = np_utils.to_categorical(y_train, nb_classes)
    y_test  = np_utils.to_categorical(y_test, nb_classes)

    # save datasets as np.ndarray class format files
    np.save('X_train', X_train)
    np.save('y_train', y_train)
    np.save('X_test' , X_test)
    np.save('y_test' , y_test)

mobilenetv2.py

MobileNetV2のアーキテクチャを定義してそのインスタンスを返すためのファイルです.
Kerasではkeras.applications.mobilenetv2.MobileNetV2で、定義ずみアーキテクチャの利用が可能なのですが, CIFAR-10, CIFAR-100の画像データは一片が32 pixelと非常に小さく、一辺が224 pixelで構成されるImageNet用に書かれている原論文のモデルでは, うまく学習ができません. そのため, Githubの実装を参考にして, アーキテクチャを作りました.

import os
import warnings
import numpy as np
from keras.layers import Input, Activation, Conv2D, Dense, Dropout, BatchNormalization, ReLU, DepthwiseConv2D, GlobalAveragePooling2D, GlobalMaxPooling2D, Add
from keras.models import Model
from keras import regularizers

# define the filter size
def _make_divisible(v, divisor, min_value=None):
    if min_value is None:
        min_value = divisor
    new_v = max(min_value, int(v + divisor / 2) // divisor * divisor)
    # Make sure that round down does not go down by more than 10%.
    if new_v < 0.9 * v:
        new_v += divisor
    return new_v


# define the calcuration of each 'Res_Block'
def _inverted_res_block(inputs, expansion, stride, alpha, filters, block_id):
    prefix = 'block_{}_'.format(block_id)

    in_channels = inputs._keras_shape[-1]
    pointwise_conv_filters = int(filters * alpha)
    pointwise_filters = _make_divisible(pointwise_conv_filters, 8)
    x = inputs

    # Expand
    if block_id:
        x = Conv2D(expansion * in_channels, kernel_size=1, strides=1, padding='same', use_bias=False, activation=None, kernel_initializer="he_normal", kernel_regularizer=regularizers.l2(4e-5), name=prefix + 'expand')(x)
        x = BatchNormalization(epsilon=1e-3, momentum=0.999, name=prefix + 'expand_BN')(x)
        x = ReLU(6., name=prefix + 'expand_relu')(x)
    else:
        prefix = 'expanded_conv_'

    # Depthwise
    x = DepthwiseConv2D(kernel_size=3, strides=stride, activation=None, use_bias=False, padding='same', kernel_initializer="he_normal", depthwise_regularizer=regularizers.l2(4e-5), name=prefix + 'depthwise')(x)
    x = BatchNormalization(epsilon=1e-3, momentum=0.999, name=prefix + 'depthwise_BN')(x)
    x = ReLU(6., name=prefix + 'depthwise_relu')(x)

    # Project
    x = Conv2D(pointwise_filters, kernel_size=1, strides=1, padding='same', use_bias=False, activation=None, kernel_initializer="he_normal", kernel_regularizer=regularizers.l2(4e-5), name=prefix + 'project')(x)
    x = BatchNormalization(epsilon=1e-3, momentum=0.999, name=prefix + 'project_BN')(x)


    if in_channels == pointwise_filters and stride == 1:
        return Add(name=prefix + 'add')([inputs, x])
    return x

# build MobileNetV2 models
def MobileNetV2(input_shape=(32, 32, 3),
                alpha=1.0,
                depth_multiplier=1,
                include_top=True,
                pooling=None,
                classes=10):

    # fileter size (first block)
    first_block_filters = _make_divisible(32 * alpha, 8)
    # input shape  (first block)
    img_input = Input(shape=input_shape)

    # model architechture
    x = Conv2D(first_block_filters, kernel_size=3, strides=1, padding='same', use_bias=False, kernel_initializer="he_normal", kernel_regularizer=regularizers.l2(4e-5), name='Conv1')(img_input)
    #x = BatchNormalization(epsilon=1e-3, momentum=0.999, name='bn_Conv1')(x)
    #x = ReLU(6., name='Conv1_relu')(x)

    x = _inverted_res_block(x, filters=16,  alpha=alpha, stride=1, expansion=1, block_id=0 )

    x = _inverted_res_block(x, filters=24,  alpha=alpha, stride=1, expansion=6, block_id=1 )
    x = _inverted_res_block(x, filters=24,  alpha=alpha, stride=1, expansion=6, block_id=2 )

    x = _inverted_res_block(x, filters=32,  alpha=alpha, stride=2, expansion=6, block_id=3 )
    x = _inverted_res_block(x, filters=32,  alpha=alpha, stride=1, expansion=6, block_id=4 )
    x = _inverted_res_block(x, filters=32,  alpha=alpha, stride=1, expansion=6, block_id=5 )

    x = _inverted_res_block(x, filters=64,  alpha=alpha, stride=2, expansion=6, block_id=6 )
    x = _inverted_res_block(x, filters=64,  alpha=alpha, stride=1, expansion=6, block_id=7 )
    x = _inverted_res_block(x, filters=64,  alpha=alpha, stride=1, expansion=6, block_id=8 )
    x = _inverted_res_block(x, filters=64,  alpha=alpha, stride=1, expansion=6, block_id=9 )

    x = _inverted_res_block(x, filters=96,  alpha=alpha, stride=1, expansion=6, block_id=10)
    x = _inverted_res_block(x, filters=96,  alpha=alpha, stride=1, expansion=6, block_id=11)
    x = _inverted_res_block(x, filters=96,  alpha=alpha, stride=1, expansion=6, block_id=12)

    x = _inverted_res_block(x, filters=160, alpha=alpha, stride=2, expansion=6, block_id=13)
    x = _inverted_res_block(x, filters=160, alpha=alpha, stride=1, expansion=6, block_id=14)
    x = _inverted_res_block(x, filters=160, alpha=alpha, stride=1, expansion=6, block_id=15)
    x = Dropout(rate=0.25)(x)

    x = _inverted_res_block(x, filters=320, alpha=alpha, stride=1, expansion=6, block_id=16)
    x = Dropout(rate=0.25)(x)

    # define fileter size (last block)
    if alpha > 1.0:
        last_block_filters = _make_divisible(1280 * alpha, 8)
    else:
        last_block_filters = 1280


    x = Conv2D(last_block_filters, kernel_size=1, use_bias=False, kernel_initializer="he_normal", kernel_regularizer=regularizers.l2(4e-5), name='Conv_1')(x)
    x = BatchNormalization(epsilon=1e-3, momentum=0.999, name='Conv_1_bn')(x)
    x = ReLU(6., name='out_relu')(x)
    
    # top layer ("use" or "not use" FC)
    if include_top:
        x = GlobalAveragePooling2D(name='global_average_pool')(x)
        x = Dense(classes, activation='softmax', use_bias=True, name='Logits')(x)
    else:
        if pooling == 'avg':
            x = GlobalAveragePooling2D()(x)
        elif pooling == 'max':
            x = GlobalMaxPooling2D()(x)

    # create model of MobileNetV2 (for CIFAR-10)
    model = Model(inputs=img_input, outputs=x, name='mobilenetv2_cifar10')
    return model

main.py

各ファイルを統合して, 学習を実行するメインファイルです. このファイルによる実行がトリガーとなって, 各処理が進みます.

import os
import h5py
import numpy as np
import keras.backend.tensorflow_backend as KTF
import tensorflow as tf
from keras.models import Sequential, Model
from keras.applications.mobilenetv2 import MobileNetV2
from keras.engine.topology import Input
from keras.layers.core import Dense, Dropout, Activation, Flatten
from keras.layers.normalization import BatchNormalization
from keras.layers import Conv2D, SeparableConv2D, MaxPooling2D
from keras.optimizers import SGD, Adam, RMSprop
from keras.preprocessing.image import ImageDataGenerator
from keras.callbacks import EarlyStopping, ModelCheckpoint, TensorBoard, CSVLogger, LearningRateScheduler
from datetime import datetime

from datafeeders.datafeeder import DataFeeder
from app.mobilenet_v2 import MobileNetV2

# set meta params
batch_size = 128
nb_classes = 10
nb_epoch   = 300
nb_data    = 32*32

# set meta info
log_dir         = '../train_log/mobilenet_v2-like5_log'
dataset_dir     = '../../CIFAR-10/datasets/dataset_norm'
model_name      = 'mobilenet_v2-like5__' + datetime.now().strftime("%Y-%m-%d-%H-%M-%S")
model_arch_path = os.path.join(log_dir, (model_name + '_arch.png'))
model_cp_path   = os.path.join(log_dir, (model_name + '_checkpoint.h5'))
model_csv_path  = os.path.join(log_dir, (model_name + '_csv.csv'))


# load data
DF = DataFeeder(dataset_dir)
X_train, y_train, X_test, y_test = DF.X_train, DF.y_train, DF.X_test, DF.y_test

# data augumatation
datagen = ImageDataGenerator(
        width_shift_range=0.2, 
        height_shift_range=0.2,
        vertical_flip=False,
        horizontal_flip=True
        )
datagen.fit(X_train)


# build model
model = MobileNetV2(input_shape=X_train.shape[1:], include_top=True, alpha=1.0)
model.summary()
print('Model Name: ', model_name)

# save model architechture plot (.png)
from keras.utils import plot_model
plot_model(model, to_file=model_arch_path, show_shapes=True)

# set learning rate
learning_rates=[]
for i in range(5):
    learning_rates.append(2e-2)
for i in range(50-5):
    learning_rates.append(1e-2)
for i in range(100-50):
    learning_rates.append(8e-3)
for i in range(150-100):
    learning_rates.append(4e-3)
for i in range(200-150):
    learning_rates.append(2e-3)
for i in range(300-200):
    learning_rates.append(1e-3)

# set callbacks
callbacks = []
callbacks.append(TensorBoard(log_dir=log_dir, histogram_freq=1))
callbacks.append(ModelCheckpoint(model_cp_path, monitor='val_loss', save_best_only=True))
callbacks.append(LearningRateScheduler(lambda epoch: float(learning_rates[epoch])))
callbacks.append(CSVLogger(model_csv_path)) 

# compile & learning model
model.compile(loss='categorical_crossentropy', optimizer=SGD(lr=2e-2, momentum=0.9, decay=0.0, nesterov=False), metrics=['accuracy'])
history = model.fit_generator(
              datagen.flow(X_train, y_train, batch_size=batch_size),
              steps_per_epoch=len(X_train) / batch_size,
              epochs=nb_epoch,
              verbose=1,
              callbacks=callbacks,
              validation_data=(X_test, y_test))
   
# validation
val_loss, val_acc = model.evaluate(X_test, y_test, verbose=0)
print('Model Name: ', model_name)
print('Test loss     : {:.5f}'.format(val_loss))
print('Test accuracy : {:.5f}'.format(val_acc))

# save model "INSTANCE"
f1_name = model_name + '_instance'
f1_path = os.path.join(log_dir, f1_name) + '.h5'
model.save(f1_path)

# save model "WEIGHTs"
f2_name = model_name + '_weights'
f2_path = os.path.join(log_dir, f2_name) + '.h5'
model.save_weights(f2_path)

# save model "ARCHITECHTURE"
f3_name = model_name + '_architechture'
f3_path = os.path.join(log_dir, f3_name) + '.json'
json_string = model.to_json()
with open(f3_path, 'w') as f:
    f.write(json_string)

# save plot figure
from utils.plot_log import save_plot_log
save_plot_log(log_dir, model_csv_path, index='acc')
save_plot_log(log_dir, model_csv_path, index='loss')

3. 学習結果

最終的な Validation Accuracyは 91% くらいでした。パラメータ数が250万弱なのでなかなかの性能。
学習率のスケジューリングによって精度が±5％動くので、難しいです。
※lossグラフの縦軸に「crossentropy」とありますが、多変量なので正しくは「categorical crossentropy」です。さーせん。

f:id:yumaloop:20180910124348p:plain — acc

f:id:yumaloop:20180910124359p:plain — loss

2018-08-15

Keras/Tensorflow : CIFAR-10のVGG-likeなアーキテクチャを作った.

Python Tensorflow 機械学習 Keras

f:id:yumaloop:20181201150912p:plain — VGG

1. 動作環境

OS: Ubuntu 16.04

Package             Version  
------------------- -------
python              3.5.0
tensorboard         1.9.0    
tensorflow          1.9.0    
h5py                2.8.0    
Keras               2.2.2    
Keras-Applications  1.0.4    
Keras-Preprocessing 1.0.2

2. プログラム

import os
import numpy as np
import keras.backend.tensorflow_backend as KTF
import tensorflow as tf
from keras.models import Sequential
from keras.layers.core import Dense, Dropout, Activation, Flatten
from keras.layers.normalization import BatchNormalization
from keras.layers import Conv2D, MaxPooling2D
from keras.optimizers import SGD, Adam
from keras.preprocessing.image import ImageDataGenerator
from keras.callbacks import EarlyStopping, ModelCheckpoint, TensorBoard, CSVLogger
from datetime import datetime

# Set Meta Parameters & Information
img_row = 32.0
img_col = 32.0
img_ch  = 255.0

batch_size = 128
nb_classes = 10
nb_epoch   = 200
nb_data    = 32*32

log_dir        = '../train_log/vgg-like_log'
dataset_dir    = '../../CIFAR-10/datasets'
model_name     = 'vgg-like__' + datetime.now().strftime("%Y-%m-%d-%H-%M-%S")
model_cp_path  = os.path.join(log_dir, (model_name + '_checkpoint.h5'))
model_csv_path = os.path.join(log_dir, (model_name + '_csv.csv'))

# Load CIFAR-10 Dataset
argmax_ch  = 255.0
(X_train, y_train), (X_test, y_test) = cifar10.load_data()

# convert pixel value range 0.0-255.0 to 0.0-1.0
X_train = X_train.astype('float32') / img_ch
X_test  = X_test.astype('float32')  / img_ch

# convert class label (0-9) to one-hot encoding format
y_train = np_utils.to_categorical(y_train, nb_classes)
y_test  = np_utils.to_categorical(y_test, nb_classes)

# save datasets as "np.ndarray" format files
np.save('X_train', X_train)
np.save('y_train', y_train)
np.save('X_test' , X_test)
np.save('y_test' , y_test)

# Data Augumatation
datagen = ImageDataGenerator(
        featurewise_center=True,
        featurewise_std_normalization=True,
        rotation_range=20,
        width_shift_range=0.2,
        height_shift_range=0.2,
        horizontal_flip=True
        vertical_flip=False)
datagen.fit(X_train)



old_session = KTF.get_session()

with tf.Graph().as_default():
    session = tf.Session('')
    KTF.set_session(session)
    KTF.set_learning_phase(1)

    # build model
    model = Sequential()
    with tf.name_scope('inference') as scope:
        model.add(Conv2D(64, (3, 3), padding='same', input_shape=X_train.shape[1:]))
        model.add(BatchNormalization())
        model.add(Activation('relu'))
        model.add(Conv2D(64, (3, 3), padding='same'))
        model.add(BatchNormalization())
        model.add(Activation('relu'))
        model.add(MaxPooling2D(pool_size=(2, 2)))
        model.add(Dropout(0.25))

        model.add(Conv2D(128, (3, 3), padding='same'))
        model.add(BatchNormalization())
        model.add(Activation('relu'))
        model.add(Conv2D(128, (3, 3), padding='same'))
        model.add(BatchNormalization())
        model.add(Activation('relu'))
        model.add(MaxPooling2D(pool_size=(2, 2)))
        model.add(Dropout(0.25))

        model.add(Conv2D(256, (3, 3), padding='same'))
        model.add(BatchNormalization())
        model.add(Activation('relu'))
        model.add(Conv2D(256, (3, 3), padding='same'))
        model.add(BatchNormalization())
        model.add(Activation('relu'))
        model.add(Conv2D(256, (3, 3), padding='same'))
        model.add(BatchNormalization())
        model.add(Activation('relu'))
        model.add(Conv2D(256, (3, 3), padding='same'))
        model.add(BatchNormalization())
        model.add(Activation('relu'))
        model.add(MaxPooling2D(pool_size=(2, 2)))
        model.add(Dropout(0.25))

        model.add(Flatten())
        model.add(Dense(1024))
        model.add(Activation('relu'))
        model.add(Dropout(0.5))
        model.add(Dense(1024))
        model.add(Activation('relu'))
        model.add(Dropout(0.5))
        model.add(Dense(nb_classes))
        model.add(Activation('softmax'))
    model.summary()

    # set callbacks
    cp_cb  = ModelCheckpoint(model_cp_path, monitor='val_loss', save_best_only=True)
    tb_cb  = TensorBoard(log_dir=log_dir, histogram_freq=1)
    csv_cb = CSVLogger(model_csv_path) 
    callbacks = [cp_cb, tb_cb, csv_cb]

    # compile model
    model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
    history = model.fit_generator(
        datagen.flow(X_train, y_train, batch_size=batch_size),
        steps_per_epoch=len(X_train) / batch_size,
        epochs=nb_epoch,
        verbose=1,
        callbacks=callbacks,
        validation_data=(X_test, y_test))
    
    # validation
    score = model.evaluate(X_test, y_test, verbose=0)
    print('val score    : ', score[0])
    print('val accuracy : ', score[1])


# save model "INSTANCE"
f1_name = model_name + '_instance'
f1_path = os.path.join(log_dir, f1_name) + '.h5'
model.save(f1_path)

# save model "WEIGHTs"
f2_name = model_name + '_weights'
f2_path = os.path.join(log_dir, f2_name) + '.h5'
model.save_weights(f2_path)

# save model "ARCHITECHTURE"
f3_name = model_name + '_architechture'
f3_path = os.path.join(log_dir, f3_name) + '.json'
json_string = model.to_json()
with open(f3_path, 'w') as f:
    f.write(json_string)


# end of session
KTF.set_session(old_session)

3.アーキテクチャ

Kerasのkeras.models.Model()クラスもつ属性(attribute)である"summary()"を使う.
プログラム上のmodel.summary()で、標準出力にモデルの構造(architechture)の要約情報が表示される.

Modelクラス (functional API) - Keras Documentation

_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
conv2d_1 (Conv2D)            (None, 32, 32, 64)        1792      
_________________________________________________________________
batch_normalization_1 (Batch (None, 32, 32, 64)        256       
_________________________________________________________________
activation_1 (Activation)    (None, 32, 32, 64)        0         
_________________________________________________________________
conv2d_2 (Conv2D)            (None, 32, 32, 64)        36928     
_________________________________________________________________
batch_normalization_2 (Batch (None, 32, 32, 64)        256       
_________________________________________________________________
activation_2 (Activation)    (None, 32, 32, 64)        0         
_________________________________________________________________
max_pooling2d_1 (MaxPooling2 (None, 16, 16, 64)        0         
_________________________________________________________________
dropout_1 (Dropout)          (None, 16, 16, 64)        0         
_________________________________________________________________
conv2d_3 (Conv2D)            (None, 16, 16, 128)       73856     
_________________________________________________________________
batch_normalization_3 (Batch (None, 16, 16, 128)       512       
_________________________________________________________________
activation_3 (Activation)    (None, 16, 16, 128)       0         
_________________________________________________________________
conv2d_4 (Conv2D)            (None, 16, 16, 128)       147584    
_________________________________________________________________
batch_normalization_4 (Batch (None, 16, 16, 128)       512       
_________________________________________________________________
activation_4 (Activation)    (None, 16, 16, 128)       0         
_________________________________________________________________
max_pooling2d_2 (MaxPooling2 (None, 8, 8, 128)         0         
_________________________________________________________________
dropout_2 (Dropout)          (None, 8, 8, 128)         0         
_________________________________________________________________
conv2d_5 (Conv2D)            (None, 8, 8, 256)         295168    
_________________________________________________________________
batch_normalization_5 (Batch (None, 8, 8, 256)         1024      
_________________________________________________________________
activation_5 (Activation)    (None, 8, 8, 256)         0         
_________________________________________________________________
conv2d_6 (Conv2D)            (None, 8, 8, 256)         590080    
_________________________________________________________________
batch_normalization_6 (Batch (None, 8, 8, 256)         1024      
_________________________________________________________________
activation_6 (Activation)    (None, 8, 8, 256)         0         
_________________________________________________________________
conv2d_7 (Conv2D)            (None, 8, 8, 256)         590080    
_________________________________________________________________
batch_normalization_7 (Batch (None, 8, 8, 256)         1024      
_________________________________________________________________
activation_7 (Activation)    (None, 8, 8, 256)         0         
_________________________________________________________________
conv2d_8 (Conv2D)            (None, 8, 8, 256)         590080    
_________________________________________________________________
batch_normalization_8 (Batch (None, 8, 8, 256)         1024      
_________________________________________________________________
activation_8 (Activation)    (None, 8, 8, 256)         0         
_________________________________________________________________
max_pooling2d_3 (MaxPooling2 (None, 4, 4, 256)         0         
_________________________________________________________________
dropout_3 (Dropout)          (None, 4, 4, 256)         0         
_________________________________________________________________
flatten_1 (Flatten)          (None, 4096)              0         
_________________________________________________________________
dense_1 (Dense)              (None, 1024)              4195328   
_________________________________________________________________
activation_9 (Activation)    (None, 1024)              0         
_________________________________________________________________
dropout_4 (Dropout)          (None, 1024)              0         
_________________________________________________________________
dense_2 (Dense)              (None, 1024)              1049600   
_________________________________________________________________
activation_10 (Activation)   (None, 1024)              0         
_________________________________________________________________
dropout_5 (Dropout)          (None, 1024)              0         
_________________________________________________________________
dense_3 (Dense)              (None, 10)                10250     
_________________________________________________________________
activation_11 (Activation)   (None, 10)                0         
=================================================================
Total params: 7,586,378
Trainable params: 7,583,562
Non-trainable params: 2,816
_________________________________________________________________

4.学習結果

validation accuracy は 90% くらい

f:id:yumaloop:20180822170332p:plain

f:id:yumaloop:20180822170409p:plain

1. KL-divergenceとは?

1.1 定義

1.2 基本的な性質

1.3 KL-divergenceは距離なのか？

2. 諸量との関係

2.1 KL-divergenceと相互情報量

2.2 KL-divergenceと対数尤度比

2.3 KL-divergenceとFisher情報量

3. 参考書籍

1. 推定のモチベーション

2. 最尤推定（ML）

2.1 パラメータ の対数尤度

2.2 パラメータの推定（最尤推定量）

3. 事後確率最大化推定（MAP）

3.1 パラメータ の事後確率

3.2 パラメータの推定（MAP推定量）

4. ベイズ推定（Bayse）

4.1 パラメータ の平均対数尤度

4.2 確率分布 の推定（予測分布）

5. 参考図書

1. 確率収束 : convergence in probability

2. p次平均収束 : convergence in mean

3. 分布収束（法則収束）: convergence in distribution (law)

4. 概収束 : almost sure convergence

参考文献

ベイズ推定とは？

1. 事後分布

2. 予測分布

3. 損失関数

4. 誤差関数

5. 汎化損失・経験損失

1. 推測統計学の概念・用語

1.1. 推測統計学とそのモチベーション

1.2. 概念・用語の定義

2. 統計的仮説検定

2.1. 母数と標本

2.2. 仮説と検定

2.3. 有意水準 , P値

3. 正規母集団についての仮説検定

3.1. 1標本の正規母集団

3.1.1. 母平均 の検定（母分散 が既知）

3.1.2. 母平均 の検定（母分散 が未知）

3.1.3. 母分散 の検定

3.2 2標本の正規母集団（母分散が等しい）

3.2.1 母平均 の同等性の検定（母分散 が既知）

3.2.2 母平均 の同等性の検定（母分散 が未知） "t検定"

3.3 2標本の正規母集団（母分散が等しくない）

3.3.1 母平均 の同等性の検定（母分散 が既知）

3.3.2 母平均 の同等性の検定（母分散 が未知）

3.3.3 母分散 の同等性の検定 "F検定"

参考

0. 前提

1. 動作環境

2. プログラム

import_CIFAR-10.py

mobilenetv2.py

main.py

3. 学習結果

1. 動作環境

2. プログラム

3.アーキテクチャ

4.学習結果

2.1 パラメータ $w$ の対数尤度 $\ell(w|\cdot)$

3.1 パラメータ $w$ の事後確率 $p(w|\cdot)$

4.1 パラメータ $w$ の平均対数尤度 $K$

4.2 確率分布 $q(x)$ の推定（予測分布）

1. 事後分布 $p(w|X_n)$

2. 予測分布 $p^{*}(x)$

3. 損失関数 $L(~)$

4. 誤差関数 $K(~)$

1.1.　推測統計学とそのモチベーション

1.2.　概念・用語の定義

2.1.　母数と標本

2.2.　仮説と検定

2.3.　有意水準 $\alpha$ , P値

3.1.　1標本の正規母集団

3.1.1.　母平均 $\mu$ の検定（母分散 ${\sigma}^{2}$ が既知）

3.1.2.　母平均 $\mu$ の検定（母分散 ${\sigma}^{2}$ が未知）

3.1.3.　母分散 ${\sigma}^2$ の検定

3.2　2標本の正規母集団（母分散が等しい）

3.2.1　母平均 $\mu_1,~\mu_2$ の同等性の検定（母分散 ${{\sigma}^2}$ が既知）

3.2.2　母平均 $\mu_1,~\mu_2$ の同等性の検定（母分散 ${{\sigma}^2}$ が未知）　"t検定"

3.3　2標本の正規母集団（母分散が等しくない）

3.3.1　母平均 $\mu_1,~\mu_2$ の同等性の検定（母分散 ${{\sigma_1}^2},~{{\sigma_2}^2}$ が既知）

3.3.2　母平均 $\mu_1,~\mu_2$ の同等性の検定（母分散 ${{\sigma_1}^2},~{{\sigma_2}^2}$ が未知）

3.3.3　母分散 ${\sigma_1}^2,~{\sigma_2}^2$ の同等性の検定　"F検定"