正規分布の推定

概要

平均 \(\mu\)、分散 \(\sigma^2\) をパラメータとする正規分布 \({\rm Norm}(\mu,\sigma^2)\) にもとづいて確率変数 \(x\) が観測されるとき、その確率密度関数は以下のように表される。\[ P(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} {\rm exp}\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\} \]

尤度関数 =正規分布の同時確率

この確率分布にもとづいて \(n\) 個のデータ \((x_1, x_2, \ldots, x_n) = \vector{x}\) の組み合わせが観測される同時確率はそれぞれの \(x_i\) が独立であれば以下のように表される。\[ \begin{eqnarray} P(\vector{x}| \mu,\sigma^2) & = & \prod_{i=1}^n P(x_i|\mu,\sigma^2) \nonumber \\ & = & \left(\frac{1}{\sqrt{2 \pi \sigma^2}}\right)^n \exp \left\{-\frac{1}{2\sigma^2} \sum_{i=1}^n (x_i-\mu)^2 \right\} \label{likelihood_function} \end{eqnarray} \]

実際にデータ \((x_1, x_2, \ldots, x_n) = \vector{x}\) を観測したとき \(P(x_1,\cdots,x_n|\mu,\sigma^2)\) は \(\vector{x}\) に対する \(\mu\) の確度を表す尤度関数となる。ここで対数尤度関数は以下のように表される。\[ \begin{equation} \label{log_likelihood} \log P(\vector{x}| \mu,\sigma^2) = n \log \frac{1}{\sqrt{2\pi\sigma^2}} - \frac{1}{2\sigma^2} \sum_{i=1}^n (x_i-\mu)^2 \end{equation} \]

最尤推定

最尤推定は (対数) 尤度関数が最大値を取る \(\hat{\mu}_{\rm mle}\) を求める。\(\hat{\mu}_{\rm mle}\) は対数尤度関数 (\(\ref{log_likelihood}\)) が極値となる位置の \(\mu\) である。したがって対数尤度関数の微分値が 0 となる位置の \(\mu\) を求める。\[ \frac{\partial}{\partial \mu} \left\{ n \log \frac{1}{\sqrt{2\pi\sigma^2}} - \sum_{i=1}^n \frac{(x_i-\mu)^2}{2\sigma^2} \right\} = \frac{n\mu}{\sigma^2} - \sum_{i=0}^n \frac{x_i}{\sigma^2} = 0 \] したがって \[ \hat{\mu}_{\rm mle} = \frac{1}{n}\sum_{i=1}^n x_i \]

正規分布に対して最尤推定で求めた \(\hat{\mu}_{\rm mle}\) は相加平均と等価である。

ベイズ推定

正規分布は平均 \(\mu\) と分散 \(\sigma^2\) の 2 パラメータによって決定する。ベイズの定理で事後確率を求めるにはそれらのどちらか一つに対してしか行うことができない。ここでは平均 \(\mu\) の分布と分散 \(\sigma^2\) の分布とを個別に求める。

正規分布に従うある標本を観測する場合、事前確率はそこまでに観測した標本によって得られた正規分布であり、新しく観測した標本によってそれを更新する

平均 \(\mu\) が未知の事後確率

分散 \(\sigma^2\) が分かっているが平均 \(\mu\) が未知の正規分布 \({\rm Norm}(\mu;\sigma^2)\) に従う標本 \(\vector{x} = (x_1,\cdots,x_n)\) を観測した場合について考える。事後確率を導く上で必要な数式変換について先に述べておく。

平均との差の自乗和

\[ \begin{equation} \sum_{i=1}^n (x_i - \mu)^2 = \sum_{i=1}^n (x_i - \bar{x})^2 + n(\bar{x} - \mu)^2 \label{sum_of_differences_from_the_mean} \end{equation} \] ここで標本平均 \(\displaystyle \bar{x} = \frac{1}{n} \sum_{i=1}^n x_i\) とする。

独立かつ同一分布に従う確率変数の集合 \(\vector{x} = (x_1, \cdots, x_n)\) が分散 \(\sigma^2\) が既知の正規分布 \(\vector{x} \sim {\rm Norm}(\mu,\sigma^2)\) に従う場合、共役事前分布も正規分布である。

ここで数式の簡略化のため \(\tau = \frac{1}{\sigma^2}\) とする (分散の逆数である \(\tau\) は精度 (precision) と呼ばれる)。尤度関数は正規分布の \(n\) 個の同時確率として式 (\(\ref{mu_likelihood}\)) のように表される。\[ \begin{eqnarray} P(\vector{x}|\mu,\tau) & = & \prod_{i=1}^n \sqrt{\frac{\tau}{2\pi}} {\rm exp}\left\{-\frac{\tau}{2}(x_i-\mu)^2\right\} \nonumber \\ & = & \left(\frac{\tau}{2\pi}\right)^{\frac{n}{2}} {\rm exp}\left\{-\frac{\tau}{2}\sum_{i=1}^n (x_i-\mu)^2\right\} \nonumber \\ & = & \left(\frac{\tau}{2\pi}\right)^{\frac{n}{2}} {\rm exp}\left[-\frac{\tau}{2}\left\{\sum_{i=1}^n (x_i-\bar{x})^2 + n(\bar{x}-\mu)^2\right\}\right] \label{mu_likelihood} \end{eqnarray} \]

式 (\(\ref{mu_likelihood}\)) より事後確率は以下のように表される。\[ \begin{eqnarray} P(\mu|\vector{x}) & \propto & P(\vector{x}|\mu) P(\mu) \nonumber \\ & = & \left(\frac{\tau}{2\pi}\right)^{\frac{n}{2}} {\rm exp}\left[-\frac{\tau}{2}\left\{\sum_{i=1}^n (x_i-\bar{x})^2 + n(\bar{x}-\mu)^2\right\}\right] \sqrt{\frac{\tau_0}{2\pi}} {\rm exp}\left\{-\frac{\tau_0}{2}(\mu-\mu_0)^2\right\} \nonumber \\ & \propto & {\rm exp}\left[-\frac{1}{2}\left\{\tau\left(\sum_{i=1}^n(x_i-\bar{x})^2+n(\bar{x}-\mu)^2\right)+\tau_0(\mu-\mu_0)^2\right\}\right] \nonumber \\ & \propto & {\rm exp}\left[-\frac{1}{2}\left\{n\tau(\bar{x}-\mu)^2+\tau_0(\mu-\mu_0)^2\right\}\right] \nonumber \\ & = & {\rm exp}\left\{-\frac{1}{2}(n\tau+\tau_0)\left(\mu-\frac{n\tau\bar{x}+\tau_0\mu_0}{n\tau+\tau_0}\right)^2+\frac{n\tau\tau_0}{n\tau+\tau_0}(\bar{x}-\mu_0)^2\right\} \nonumber \\ & \propto & {\rm exp}\left\{-\frac{1}{2}(n\tau+\tau_0)\left(\mu-\frac{n\tau\bar{x}+\tau_0\mu_0}{n\tau+\tau_0}\right)^2\right\} \label{mu_posterior} \end{eqnarray} \] 上記の変形に関して前述の二次方程式の解を使用している。また全ての定数は \(\mu\) に関与してないと想定している。結果として事後確率 (\(\ref{mu_posterior}\)) は平均 \(\frac{n\tau\bar{x}+\tau_0\mu_0}{n\tau+\tau_0}\)、精度 \(n\tau+tau_0\) の正規分布に従う。\[ P(\mu|\vector{x}) = {\rm Norm}\left(\frac{n\tau\bar{x}+\tau_0\mu_0}{n\tau+\tau_0}, \frac{1}{n\tau+\tau_0}\right) \]

分散 \(\sigma^2\) が未知の事後分布

独立かつ同一分布に従う確率変数の集合 \(\vector{x} = (x_1, \cdots, x_n)\) が平均 \(\mu\) が既知の正規分布 \(\vector{x} \sim {\rm Norm}(\mu,\sigma^2)\) に従う場合、共役事前分布は逆ガンマ関数 (inverse gamma distribution) またはScaled 逆カイ二乗分布 (scaled chi-squared distribution) となる。逆ガンマ分布を使用するのが一般的であるがここでは簡便な Scaled 逆カイ二乗分布を使用する。

Scaled 逆カイ二乗分布は以下の確率密度関数で表される。\[ f(x;\nu,\tau^2) = \frac{(\frac{\tau^2\nu}{2})^\frac{\nu}{2}}{\Gamma(\frac{\nu}{2})} \frac{{\rm exp}(-\frac{\nu\tau^2}{2x})}{x^{1+\frac{\nu}{2}}} \] これより \(\sigma^2\) に対する事前確率は以下のように表される。\[ \begin{eqnarray} P(\sigma^2|\nu_0,\sigma_0^2) & = & \frac{(\frac{\sigma_0^2\nu_0}{2})^\frac{\nu_0}{2}}{\Gamma(\frac{\nu_0}{2})} \frac{{\rm exp}(-\frac{\nu_0\sigma_0^2}{2\sigma^2})}{(\sigma^2)^{1+\frac{\nu_0}{2}}} \nonumber \\ & \propto & \frac{{\rm exp}\left(-\frac{\nu_0\sigma_0^2}{2\sigma^2}\right)}{(\sigma^2)^{1+\frac{\nu_0}{2}}} \label{sigma_prior} \end{eqnarray} \] 尤度関数 (\(\ref{likelihood_function}\)) から分散についての項をまとめると \[ \begin{eqnarray} P(\vector{x}|\mu,\sigma^2) & = & \left(\frac{1}{2\pi\sigma^2}\right)^{\frac{n}{2}} {\rm exp}\left\{-\frac{1}{2\sigma^2}\sum_{i=1}^n (x_i-\mu)^2\right\} \nonumber \\ & = & \left(\frac{1}{2\pi\sigma^2}\right)^{\frac{n}{2}} {\rm exp} \left(-\frac{S}{2\sigma^2}\right) \label{sigma_likelihood} \end{eqnarray} \] ここで \(\displaystyle S = \sum_{i=1}^n (x_i-\mu)^2\) とする。

尤度関数 (\(\ref{sigma_likelihood}\)) と事前分布 (\(\ref{sigma_prior}\)) より事後分布を求めると: \[ \begin{eqnarray} P(\sigma^2|\vector{x}) & \propto & P(\vector{x}|\sigma^2) P(\sigma^2) \nonumber \\ & \propto & \left(\frac{1}{2\pi\sigma^2}\right)^{\frac{n}{2}} {\rm exp} \left(-\frac{S}{2\sigma^2}\right) \frac{{\rm exp}\left(-\frac{\nu_0\sigma_0^2}{2\sigma^2}\right)}{(\sigma^2)^{1+\frac{\nu_0}{2}}} \nonumber \\ & = & \left(\frac{1}{2\pi\sigma^2}\right)^{\frac{n}{2}} \frac{1}{(\sigma^2)^{1+\frac{\nu_0}{2}}} {\rm exp}\left(-\frac{S + \nu_0\sigma_0^2}{2\sigma^2}\right) \nonumber \\ & = & \frac{1}{(\sigma^2)^{1+\frac{\nu_0+n}{2}}} {\rm exp}\left(-\frac{\nu_0\sigma_0^2 + S}{2\sigma^2}\right) \label{sigma_posterior} \end{eqnarray} \] 式 (\(\ref{sigma_posterior}\)) は Scaled 逆カイ二乗分布である。ここで: \[ \left\{\begin{array}{rcl} \nu_0' & = & \nu_0 + n \\ \nu_0'\sigma'{}_0^2 & = & \nu_0 \sigma_0^2 + \sum_{i=1}^n (x_i - \mu)^2 \end{array}\right. \]

共役事前分布に逆ガンマ分布を使用した場合の事後確率のパラメータは以下のようになる。\[ \left\{\begin{array}{rcl} \alpha' & = & \alpha + \frac{n}{2} \\ \beta' & = & \beta + \frac{\sum_{i=1}^n (x_i - \mu)^2}{2} \end{array}\right. \]

平均 \(\mu\) と分散 \(\sigma^2\) が未知の事後分布

数式が凄まじいので意味がわかるベイズ統計学参照。

MAP推定

MAP 推定はベイズ定理における事後確率 \(P(p|\vector{x})\) を最大化する \(\hat{\mu}_{\rm map}\) を求める。\[ P(p|\vector{x}) = \frac{P(\vector{x}|p) P(p)}{P(\vector{x})} \]

エビデンスの \(P(\vector{x})\) はすでに観測済みの \(\vector{x}\) に対する確率であり、事後確率が最大となる \(p\) の結果には関与しないため無視する。\[ P(p|\vector{x}) \propto P(\vector{x}|p) P(p) \]

ここで尤度関数 \(P(\vector{x}|p)\) は (\(\ref{likelihood_function}\)) で表される。正規表現の事前確率は平均 \(\mu\) を確率変数とするか分散 \(\sigma^2\) を確率変数とするかで変わってくる。

事前確率 = 逆ガンマ分布

事前確率 \(P(p)\) を逆ガンマ分布 (inverse gamma distribution) に従う確率分布とする。逆ガンマ分布は以下のように表される。\[ P(p; \alpha,\beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} p^{-\alpha-1} e^{-\frac{\beta}{p}}, \ \ 0 \lt x, 0 \lt \alpha, 0 \lt \beta \]

事前確率 = 正規分布

事前確率 \(P(p)\) を正規分布 \({\rm Norm}(\mu_\mu=0,\sigma_\mu^2)\) に従う確率分布とする。\[ \begin{align*} \log P(p) & = \log {\rm Norm}(x=p; \mu_\mu=0,\sigma_\mu^2) \\ & = \log \left\{\frac{1}{\sqrt{2\pi\sigma_\mu^2}} {\rm exp} \left(-\frac{p^2}{2\sigma_\mu^2}\right) \right\} \\ & = \log \frac{1}{\sqrt{2\pi\sigma_\mu^2}} - \frac{p^2}{2\sigma_\mu^2} \end{align*} \] 分布の対数をとっても事後確率が最大化する位置の \(p\) には影響しないことから計算の単純化のため対数化する。

\(p\) で偏微分して 0 となる極値を求めると MAP 推定の観点で最も確度の高い \(\hat{\mu}_{\rm map}\) が得られる。\[ \frac{\partial}{\partial p} \left\{ \sum_{i=1}^n \frac{(x_i-p)^2}{\sigma^2} + \frac{p^2}{\sigma_\mu^2} \right\} = - \frac{2}{\sigma^2}\left(\sum_{i=1}^n x_i - np\right) + \frac{2p}{\sigma_\mu^2} = 0 \\ \]

したがって \[ \begin{equation} \label{map_mu} \hat{\mu}_{\rm map} = \frac{\sigma_\mu^2}{n \sigma_\mu^2 + \sigma^2} \sum_{i=1}^n x_i \end{equation} \]

(\(\ref{map_mu}\)) より、\(n \to \infty\) の極値で \(\hat{\mu}_{\rm map} = \hat{\mu}_{\rm mle}\) となる。\[ \lim_{n \to \infty} \hat{\mu}_{\rm map} = \hat{\mu}_{\rm mle} \]

事前確率 \(P(p)\) の分散 \(\sigma_\mu^2\) を特定できない場合、\(n\) が大きくなるにつれ分散 \(\sigma_\mu^2\) と \(\sigma^2\) が連動して小さくなると想定できるため \(\sigma_\mu^2 = t\sigma^2\) とおく。\[ \hat{\mu}_{{\rm map}} \mid_{\sigma_\mu^2 = a\sigma^2} = \frac{t}{t n + 1} \sum_{i=1}^n x_i \] このとき \(t\) は \(\hat{\mu}_{\rm map}\) が \(\hat{\mu}_{\rm mle}\) に収束する速度を調整するためのパラメータとして利用することができる。\(t\) を大きくすると収束が速くなることから \(\hat{\mu}_{\rm mle}\) との乖離を小さくしたいときは \(t\) を大きく設定すると良い。

また (\(\ref{map_mu}\)) は \(\sigma_\mu^2 \to \infty\) の極値を取っても (つまり \(P(p)\) を一様分布に近づけても) \(\hat{\mu}_{\rm mle}\) と等しくなる。\[ \lim_{\sigma_\mu^2 \to \infty} \hat{\mu}_{{\rm map}} = \frac{1}{n} \sum_{i=1}^n x_i = \hat{\mu}_{\rm mle} \] つまり最尤推定は MAP 推定における \(P(x|p) P(p)\) の事前確率 \(P(p)\) が無情報事前分布 (一様分布) のケースと考えることができる。

参考文献

1. 一石賢 (2016), "意味が分かるベイズ統計学", ベレ出版 P167