ぺ ん ぎ ん の 閃 き

閃き- blog

きらびやかに、美しく、痛烈に.

指数分布とワイブル分布をPythonでプロットしてみる


 生存時間解析など、応用範囲の広い指数分布についてまとめます。指数型分布族の仲間としては、ワイブル分布・ガンマ分布の他にも、ポアソン分布、レイリー分布、ラプラス分布(両側指数分布)などがあります。

1. 前提

1.1 確率分布の定義

 指数分布確率密度関数  f と累積分布関数  F は、次のように定義されます。

  • 指数分布

$$ X \sim Exp(\lambda), ~~~ \lambda \geq 0 $$

$$ \begin{eqnarray} f(x | \lambda) & = & \lambda ~ exp \left( - \lambda x \right) ~~~ (x \geq 0) \\ F(x | \lambda) & = & 1 - exp \left( - \lambda x \right) ~~~ (x \geq 0) \\ \mathbb{E}\left[ X \right] & = & \frac{1}{\lambda} \\ \mathbb{V}\left[ X \right] & = & \frac{1}{\lambda^{2}} \end{eqnarray} $$

 また、ワイブル分布は、次のように定義されます。なお、 \Gamma(~) は、ガンマ関数*1(階乗を一般化したもの)です。

  • ワイブル分布

$$ X \sim Weibull(\alpha, \beta), ~~~ \alpha \geq 0, ~~~ \beta \geq 0 $$

$$ \begin{eqnarray} f(x | \alpha, \beta) & = & \frac{\alpha}{\beta} { \left( \frac{x}{\beta} \right) }^{\alpha - 1} \left\{ - { \left( \frac{x}{\beta} \right) }^{\alpha} \right\} ~~~ (x \geq 0) \\ F(x | \alpha, \beta) & = & 1 - exp\left\{ - {\left( \frac{x}{\beta} \right) }^{a} \right\} ~~~ (x \geq 0) \\ \mathbb{E}\left[ X \right] & = & \beta ~ \Gamma \left( 1 + \frac{1}{\alpha} \right) \\ \mathbb{V}\left[ X \right] & = & \beta^{\large \frac{2}{\alpha}} \left\{ \Gamma(1 + \frac{2}{\alpha}) - {\left( \Gamma(1 + \frac{1}{\alpha}) \right)}^{2} \right\} \end{eqnarray} $$


1. 2 指数分布とワイブル分布の関係

 ワイブル分布のパラメータをそれぞれ  \alpha=1,~ \beta=\frac{1}{\lambda}とすると指数分布に一致します。つまり、ワイブル分布は、指数分布の一般に拡張したものになります。

$$ Exp(\lambda) = Weibull(1, \frac{1}{\lambda}) $$

ちなみに、以下で定義されるガンマ分布のパラメータをそれぞれ  \alpha=1,~ \beta=\frac{1}{\lambda}とすると指数分布に一致します。

$$ Exp(\lambda) = Weibull(1, \frac{1}{\lambda}) = Gamma(1, \frac{1}{\lambda}) $$

  • ガンマ分布

$$ X \sim Gamma(\alpha, \beta) $$

$$ \begin{eqnarray} f(x | \alpha, \beta) & = & \frac{1}{\Gamma(\alpha)} \frac{1}{\beta} {\left( \frac{x}{\beta} \right) }^{\alpha - 1}~~~ (x \gt 0) \\ \mathbb{E}\left[ X \right] & = & \alpha \beta \\ \mathbb{V}\left[ X \right] & = & \alpha {\beta}^{2} \end{eqnarray} $$


2. Pythonによる実装

2.1 指数分布

python で分布を定義して、グラフをプロットしてみます。

import math
import random
import numpy as np
import matplotlib.pyplot as plt

# 指数分布

class Exp():
    def __init__(self):
        pass
    
    # 密度関数f
    def f(self, x, la=1):
        return la * math.exp(-1 * la * x)
    
    # 分布関数F
    def F(self, x, la=1):
        return 1 - math.exp(-1 * la * x)
    
    # Fの逆関数(samplingに使う)
    def invF(self, y, la=1):
        return -1 * (1 / la) * math.log(1 - y)
    
    # サンプリング
    def sampling(self, num, la):
        sample=[]
        for i in range(num):
            rand = random.random()
            sample.append(self.invF(rand, la))
        return np.array(sample)
    
    # 散布図
    def scatter_pdf(self, ax, data, la, alpha=0.4):
        pdf=[]
        for d in data:
            pdf.append(self.f(d, la))

        ax.scatter(data, pdf, alpha=alpha)
        ax.set_title("Exponential (lambda={:.2f})".format(la))
        ax.set_xlabel("value")
        ax.set_ylabel("probability density")
        ax.set_ylim(0, 1)
        
    # グラフ
    def plot_pdf(self, ax, data, la):
        pdf=[]
        for d in data:
            pdf.append(self.f(d, la))

        ax.plot(data, pdf)
        ax.set_title("Exponential (lambda={:.2f})".format(la))
        ax.set_xlabel("value")
        ax.set_ylabel("probability density")
        ax.set_ylim(0, 1)

実際にグラフを描画してみます。パラメータ  \lambda の値を変えて、分布の形状をみます。

# グラフをmatplotlibで描画してみる

data = np.linspace(0, 15, 100)

fig, ax = plt.subplots(3, 3, figsize=(18, 18))

Exp().plot_pdf(ax[0, 0], data, 0.2)
Exp().plot_pdf(ax[0, 1], data, 0.4)
Exp().plot_pdf(ax[0, 2], data, 0.6)
Exp().plot_pdf(ax[1, 0], data, 0.8)
Exp().plot_pdf(ax[1, 1], data, 1)
Exp().plot_pdf(ax[1, 2], data, 1.2)
Exp().plot_pdf(ax[2, 0], data, 1.4)
Exp().plot_pdf(ax[2, 1], data, 1.6)
Exp().plot_pdf(ax[2, 2], data, 1.8)

plt.show()
f:id:yumaloop:20190114000052p:plain:w1000

 プロットをみると、指数分布のパラメータ  \lambda は、尺度(ばらつき)を表していることがわかります。分布は  \lambda が大きくなるにつれ急峻になり、ばらつきが少なくなっていきます。平均について  \frac{1}{\lambda} \to 0 ~ (\lambda \to \infty)、分散について  \frac{1}{\lambda^{2}} \to 0 ~ (\lambda \to \infty) となることから、「指数分布はパラメータ  \lambda を上げると、平均は0に接近し、分散は限りなく小さくなる(0になる)」ということがわかります。

2.2 ワイブル分布

python で分布を定義して、グラフをプロットしてみます。

import math
import random
import numpy as np
import matplotlib.pyplot as plt

# ワイブル分布

class Weibull():
    def __init__(self):
        pass
        
    # 密度関数f
    def f(self, x, a, b):
        return (a/b) * ( (x/b)**(a - 1) ) * math.exp(-1 * (x/b)**a )
    
    # 分布関数F
    def F(self, x, a, b):
        return 1 - math.exp(-1 * (x/b)**a )
    
    # Fの逆関数
    def invF(self, y, a, b):
        return b * ( -1 * math.log(1 - y) )**(1/a)
    
    # サンプリング
    def sampling(self, num, a, b):
        sample=[]
        for i in range(num):
            rand = random.random()
            sample.append(self.invF(rand, a, b))
        return np.array(sample)
    
    # 散布図
    def scatter_pdf(self, ax, data, a, b, alpha=0.4):
        pdf=[]
        for d in data:
            pdf.append(self.f(d, a, b))

        ax.scatter(data, pdf, alpha=alpha)
        ax.set_title("Weibull (alpha={:.2f}, beta={:.2f})".format(a, b))
        ax.set_xlabel("value")
        ax.set_ylabel("probability density")
        ax.set_ylim(0, 1)
        
    # グラフ
    def plot_pdf(self, ax, data, a, b):
        pdf=[]
        for d in data:
            pdf.append(self.f(d, a, b))

        ax.plot(data, pdf)
        ax.set_title("Weibull (alpha={:.2f}, beta={:.2f})".format(a, b))
        ax.set_xlabel("value")
        ax.set_ylabel("probability density")
        ax.set_ylim(0, 1)

実際にグラフを描画してみます。パラメータ  \alpha, ~\beta の値を変えて、分布の形状をみます。

# グラフをmatplotlibで描画してみる

data = np.linspace(0, 15, 100)

fig, ax = plt.subplots(3, 3, figsize=(18, 18))

Weibull().plot_pdf(ax[0, 0], data, 1, 1)
Weibull().plot_pdf(ax[0, 1], data, 1, 2)
Weibull().plot_pdf(ax[0, 2], data, 1, 3)
Weibull().plot_pdf(ax[1, 0], data, 2, 1)
Weibull().plot_pdf(ax[1, 1], data, 2, 2)
Weibull().plot_pdf(ax[1, 2], data, 2, 3)
Weibull().plot_pdf(ax[2, 0], data, 3, 1)
Weibull().plot_pdf(ax[2, 1], data, 3, 2)
Weibull().plot_pdf(ax[2, 2], data, 3, 3)

plt.show()
f:id:yumaloop:20190114000114p:plain:w1000

 プロットをみると、ワイブル分布のパラメータ  \alpha は分布の形状を表し、パラメータ  \beta は、尺度(ばらつき)を表していることがわかります。実際、分布の形状は  \alpha によって支配されており、ワイブル分布は  \alpha=1 のとき指数分布、  \alpha=2 のときレイリー分布と一致します。分布は  \beta が大きくなるにつれ急峻になり、ばらつきが少なくなっていきます。

*1:ガンマ関数は、$$ \Gamma(n):=\int_{0}^{\infty} y^{n - 1} e^{-y} dy $$と定義され、 n自然数ならば \Gamma(n)=(n-1)! が成り立つ。

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


1. KL-divergenceとは?

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

1.1 定義

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

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


1.2 基本的な性質

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

  • 非負性)非負の値域をもつ。

$$ 0 \leq D_{KL}( Q \mid\mid P ) \leq \infty \tag{2.1} $$

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

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

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

$$ D_{KL}( Q \mid\mid P ) \neq D_{KL}( P \mid\mid Q ) \tag{2.3} $$

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

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


 例として、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(~) の公理

  1. 非負性     d(x, ~ y) \geq 0
  2. 完備性     d(x, ~ y) = 0 ~~ \Leftrightarrow ~~ x = y
  3. 対称性     d(x, ~ y) = d(y, ~ x)
  4. 三角不等式   d(x, ~ y) + d(y, ~ z) \geq d(x, ~ z)

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

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

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

f:id:yumaloop:20190115092949p:plain


2. 諸量との関係

2.1 KL-divergenceと相互情報量

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

$$ \begin{eqnarray} 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{eqnarray} $$

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

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

$$ \begin{eqnarray} 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{eqnarray} $$

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

$$ \begin{eqnarray} 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{eqnarray} $$

(例)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) の組み合わせ)推定する」というモチベーションが起こります。そのため、2つの分布のズレを測りたいとき、あるいは推定誤差を損失関数やリスク関数に組み込んでパラメータに対する最適化問題をときたいときに、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 で平均操作をおこなった値」です。

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

$$ \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 ) を対数尤度関数とする。)

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

$$ \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) を仮定したとき、パラメータ  \thetaFisher情報量  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 の近傍でもつ局所的な情報量を表していることがわかります。

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


3. 参考書籍

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

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

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

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

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

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



*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}) $$

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


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} $$


https://www.researchgate.net/profile/Eric_Van_de_Weg/publication/226192090/figure/download/fig1/AS:393798262771713@1470900230184/Probability-distributions-representing-the-principle-of-Bayesian-analysis-The-posterior.png

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

5. 参考図書

ベイズ統計の理論と方法

ベイズ統計の理論と方法

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

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

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

https://upload.wikimedia.org/wikipedia/commons/4/49/Notions-of-convergence-in-probability-theory.jpg

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\} は、確率変数  Up次平均収束する」という。  
 

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} $$ と表記する。          

参考

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

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

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

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

入門確率過程

入門確率過程

http://kristen-ressurs.no/images/Fysikk/Hilbert-space.jpg http://radhakrishna.typepad.com/.a/6a00d83453b94569e2014e87c56ff9970d-pi

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

ベイズ推定とは?

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

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

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

ベイズ統計の理論と方法

ベイズ統計の理論と方法

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

ベイズ推定の流れ」

  • 真の分布    :   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} $$

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

  1. 確率モデル  p(x|w)事前分布  \varphi(w) を設定する。
  2. 確率モデル  p(x|w) と事前分布  \varphi(w) によって、事後分布  p(w|X_{n})を導出する。
  3. 確率モデル  p(x|w) と事後分布  p(w|X_{n}) によって、予測分布  p^{*}(x) を導出する。
  4. 予測分布  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) は、その期待値ですから、下図右で青く塗りつぶされた部分の「確率重みつき面積(期待値)」となります。

https://www.science-emergence.com/media/images/src/381.png

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

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

仮説検定は、推測統計学の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 支持しない」の2者択一で判定するための方法論です。しかるに、有意水準  \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つの仮説(: 例えば「母集団の平均は1である」)に対して、「支持する or 支持しない」の判断を下す方法として、複数の検定方式が考えられます。


2.3. 有意水準  \alpha , P値

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

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\}

参考

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

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

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

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

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

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

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

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

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

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」の導入などが特徴的です。

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