【E資格対策】パラメータ推定手法

最尤推定

最尤推定は, 事前確率を考慮せず, 観測データが得られる確率(尤度)のみを利用して確率モデルのパラメータを点推定する手法である.

未知のデータ生成分布 Pdata(X)P_{data}(X) から独立に抽出された mm 個の観測データ X={x1,x2,,xm}X = \{x_1, x_2, \ldots, x_m\} について, その任意の要素 xix_i が確率モデル Pmodel(xi;θ)P_{model}(x_i ; \theta) によって真の確率 Pdata(xi)P_{data}(x_i) に近似されるとする. このとき, パラメータ θ\theta の最尤推定量 θ^\hat{\theta} は次式で定義される.

θ^=argmaxθPmodel(X;θ)=argmaxθi=1mPmodel(xi;θ)\hat{\theta} = \underset{\theta}{\mathrm{argmax}} \, P_{model}(X ; \theta) = \underset{\theta}{\mathrm{argmax}} \prod_{i=1}^{m} P_{model}(x_i ; \theta)

右辺の Pmodel(xi;θ)P_{model}(x_i ; \theta) は尤度関数と呼ばれる.

ここで問題なのは, 右辺が確率の総積になっていることである. 確率を何度も掛け合わせると非常に小さな値になってしまい, 数値計算上の不安定性が生じる可能性がある.

そこで, 尤度関数の対数を取った対数尤度関数を用いることで, 積を和に変換する :

θ^=argmaxθi=1mlogPmodel(xi;θ)\hat{\theta} = \underset{\theta}{\mathrm{argmax}} \sum_{i=1}^{m} \log P_{model}(x_i ; \theta)

この式全体を mm で割れば, 訓練データによって定義される経験分布 P^data\hat{P}_{data} の期待値で表すことができる :

θ^=argmaxθEXP^data[logPmodel(x;θ)](1)\hat{\theta} = \underset{\theta}{\mathrm{argmax}} \, \mathbb{E}_{X \sim \hat{P}_{data}} [\log P_{model}(x ; \theta)] \tag{1}

最尤推定を解釈する一つの方法に, 「最尤推定は訓練集合で定義される経験分布 P^data\hat{P}_{data} とモデル分布の間の差を最小化するとみなす」という方法がある. このとき, 二つの分布間の差をKLダイバージェンスで測定する :

DKL(P^dataPmodel)=xP^data(x)logP^data(x)Pmodel(x)=xP^data(x)(logP^data(x)logPmodel(x))=EXP^data[logP^data(x)](データ生成分布のエントロピー)EXP^data[logPmodel(x)]対数尤度の期待値\begin{align*} D_{KL}(\hat{P}_{data} || P_{model}) & = \sum_{x} \hat{P}_{data}(x) \log \frac{\hat{P}_{data}(x)}{P_{model}(x)} \\ &= \sum_{x} \hat{P}_{data}(x) \left( \log \hat{P}_{data}(x) - \log P_{model}(x) \right) \\ &= \underbrace{\mathbb{E}_{X \sim \hat{P}_{data}} [\log \hat{P}_{data}(x)]}_{- \text{(データ生成分布のエントロピー)}} - \underbrace{\mathbb{E}_{X \sim \hat{P}_{data}} [\log P_{model}(x)]}_{\text{対数尤度の期待値}} \end{align*}

ここで, 第1項の E[logP^data(x)]\mathbb{E} \left[ \log \hat{P}_{data}(x) \right] はデータ生成過程のみに依存する定数であり, 最適化すべきモデルパラメータ θ\theta には依存しない. したがって, KLダイバージェンス DKL(P^dataPmodel)D_{KL}(\hat{P}_{data} || P_{model}) を最小化することは, 第2項の E[logPmodel(x)]-\mathbb{E} \left[ \log P_{model}(x) \right] を最小化すること, すなわち E[logPmodel(x)]\mathbb{E} \left[ \log P_{model}(x) \right] を最大化することと等価である. これは, 式 (1)(1) の対数尤度の最大化と同じことを意味する.


最尤推定では, 与えられた xx から yy を予測するために条件付き確率 P(yx;θ)P(y | x ; \theta) を推定する. ここで, XX が全ての入力, YY が観測された全ての目的関数を表す場合, 条件付き最尤推定量は以下の式で表される :

θ^=argmaxθP(YX;θ)\hat{\theta} = \underset{\theta}{\mathrm{argmax}} P(Y | X ; \theta)

ここでの事例が独立で同一の分布に従うと仮定すると,

θ^=argmaxθi=1mlogP(yixi;θ)\hat{\theta} = \underset{\theta}{\mathrm{argmax}} \sum_{i=1}^{m} \log P(y_i | x_i ; \theta)

となる.

ここで, 最尤推定における線形回帰について, モデルが単一の予測 y^\hat{y} を生成する代わりに, 条件付き確率 P(yx)P(y | x) を生成するものとする. 線形回帰アルゴリズムの導出のために, P(yx)P(y | x) が正規分布 N(y;y^(x;w),σ2)\mathcal{N}(y ; \hat{y}(x ; w), \sigma^2) に従うと定義する. すなわち, 入力 xx に対する出力 yy の確率密度関数は以下のように表される :

P(yx;w)=12πσ2exp((yy^(x;w))22σ2)P(y | x ; w) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{(y - \hat{y}(x ; w))^2}{2 \sigma^2} \right)

関数 y^(x;w)=wTx\hat{y}(x ; w) = w^T x はガウス分布の平均の予測であり, 固定値 σ2\sigma^2 は分散であると仮定する. ここで, mm 個の訓練事例 (xi,yi)(x_i, y_i) が独立にこの分布に従うとすると, その条件付き対数尤度は以下のように計算できる :

i=1mlogP(yixi;w)=i=1mlog(12πσ2exp((yiy^(xi;w))22σ2))=i=1m(log12πσ2(yiy^(xi;w))22σ2)=mlogσm2log(2π)i=1myiy^(xi;w)22σ2\begin{align*} \sum_{i=1}^{m} \log P(y_i | x_i ; w) &= \sum_{i=1}^{m} \log \left( \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{(y_i - \hat{y}(x_i ; w))^2}{2 \sigma^2} \right) \right) \\ &= \sum_{i=1}^{m} \left( \log \frac{1}{\sqrt{2 \pi \sigma^2}} - \frac{(y_i - \hat{y}(x_i ; w))^2}{2 \sigma^2} \right) \\ &= -m \log \sigma - \frac{m}{2} \log (2 \pi) - \sum_{i=1}^{m} \frac{\| y_i - \hat{y}(x_i ; w) \|^2}{2 \sigma^2} \tag{2} \end{align*}

ただし, y^i\hat{y}_iii 番目の入力 xix_i に対する線形回帰出力である.

一方, 最小二乗誤差(MSE)の式は以下のように表される :

MSE=1mi=1myiy^(xi;w)2(3)\mathrm{MSE} = \dfrac{1}{m} \sum_{i=1}^{m} \| y_i - \hat{y}(x_i ; w) \|^2 \tag{3}

(2)(2) と式 (3)(3) を比較すると, ww に関して対数尤度を最大化すると, 平均二乗誤差を最小する場合と同じパラメータの推定量 w^\hat{w} を得ることが分かる.

以上により, 最尤推定の手順として最小二乗誤差を利用できることが正当化された.

MAP(Maximum A Posteriori)推定

MAP推定は, 観測データに基づく尤度と事前分布を組み合わせて, 事後確率が最大となるパラメータの値を推定する手法である. 後述するベイズ推定のアプローチを取り入れることで, 観測データのみに依存することなく, 精度が向上する可能性が高まる.

θ\theta をパラメータ, xx を観測データ, P(θ)P(\theta) をパラメータの事前分布とすると, MAP推定量 θ^\hat{\theta} は以下の式で定義される :

θ^=argmaxθP(xθ)P(θx)=argmaxθP(xθ)+logP(θ)\hat{\theta} = \underset{\theta}{\mathrm{argmax}} \, P(x | \theta) P(\theta | x) = \underset{\theta}{\mathrm{argmax}} \, P(x | \theta) + \log P(\theta)

ベイズ推定

ベイズ推定は, 事前分布と観測データから得られる情報を利用して, 未知のパラメータの事後分布を導き出し, その推定を行う手法である.

パラメータを θ\theta, データ集合を XX とすると, ベイズ推定における推定値は

P(θX)=P(Xθ)P(θ)P(X)P(\theta | X) = \dfrac{P(X | \theta) P(\theta)}{P(X)}

と表せる.

ここで, P(θ)P(\theta) は事前分布, P(θX)P(\theta | X) は事後分布と呼ばれ, 一般的には両者が同じ形式となるように, 事前分布として共役事前分布を選ぶ. これにより, 事後分布の計算が容易になる.

代表的な関係は以下の通り.

  • ベルヌーイ分布や二項分布 : ベータ分布
  • ポアソン分布 : ガンマ分布
  • 正規分布(平均未知, 分散既知) : 正規分布

練習問題

問題1

ある製造工場では製品の品質検査を行っており, 生産された各製品が不良品であるかどうかは, それぞれ独立にパラメータ θ\thetaにを持つベルヌーイ分布 Bern(θ)Bern(\theta) に従う. ここで, ベルヌーイ分布のパラメータ θ\theta における事前分布をベータ分布 Beta(α,β)Beta(\alpha, \beta) とする.

新たに生産された製品の中から 2020 個の製品を無作為に抽出して1個ずつ検査したところ, 77 個が不良品であった. 事前分布のパラメータを α=4,β=6\alpha = 4, \beta = 6 としたとき, 事後分布の期待値を求めよ.


解答 (クリックで展開)

ベイズ推定の元になるベイズ推論に関する問題である.

ベイズ推論では, 観測されたデータに基づいて, 事前分布と尤度に対してベイズ則を用いることで, 事後分布を更新する. これにより, 得られたパラメータの期待値や最頻値を求めることができる.


ベルヌーイ分布の共役事前分布はベータ分布であるため, 事後分布もベータ分布 Beta(α,β)Beta(\alpha', \beta') となる.

ここで, 製品 ii の検査結果を xix_i として, xi=1x_i = 1 は不良品, xi=0x_i = 0 は良品を表すことにすると, 事後分布のパラメータは以下のように更新される(nn は製品の個数).

α=α+i=1nxiβ=β+ni=1nxi\begin{align*} \alpha' &= \alpha + \sum_{i=1}^{n} x_i \\ \beta' &= \beta + n - \sum_{i=1}^{n} x_i \end{align*}

また, ベータ分布 Beta(α,β)Beta(\alpha, \beta) の期待値は αα+β\frac{\alpha}{\alpha + \beta} である.


したがって, 事後分布のパラメータは

α=4+7=11β=6+207=19\begin{align*} \alpha' &= 4 + 7 = 11 \\ \beta' &= 6 + 20 - 7 = 19 \end{align*}

であるから, 事後分布の期待値は

E[θx]=1111+19=0.36660.367.E[\theta | x] = \frac{11}{11 + 19} = 0.3666 \cdots \approx \underline{0.367}.

参考文献