行列ガンマ分布の規格化定数の計算方法について

多次元ガウス分布の精度行列のベイズ推論を行う際にあらわれる行列ガンマ分布：
\begin{align*}
\textrm{MGam}(\Lambda|a,B)
&= C(a,B) |\Lambda|^{a-1}\exp\{-\langle B, \Lambda\rangle\}
\end{align*}
の規格化定数
\begin{align*}
C(a,B)^{-1}
&= \int_{\Lambda\in \mathrm{Sym}_+(d)} |\Lambda|^{a-1}\exp\{-\langle B, \Lambda\rangle\}\,d\Lambda
\end{align*}
の計算方法について述べます.

計算したくなる場合: 多次元ガウス分布の精度行列のベイズ推論をしたいとき

一次元ガウス分布の精度行列のベイズ推論をする際には, 事前分布としてガンマ分布が表れました.
事前分布は正のパラメータである精度の分布であったから, 予測分布(あるいは規格化定数)を求めるために での積分を行いました.

一方, 多次元においては精度行列は一般に正定値の対称行列となり, パラメータの自由度が大きくなります.
そこで, 正定値の対称行列の値域による積分を導入し, 行列ガンマ分布を定義しましょう.

【ベイズ推論】一次元ガウス分布の平均・分散についてのベイズ推論の意味

前回記事では, 一次元ガウス分布を観測モデルに持つようなベイズ推論を行うために, 事前分布はパラメータ により,
\begin{align*}
p(\lambda|v )
&= z(v)\exp\left\{v_1\ln{\lambda} -v_2\lambda +v_3\lambda\mu -\dfrac{v_4}{2}\lambda\mu^2 \right\},\\
z(v)^{-1}
&=\sqrt{\dfrac{2\pi}{v_4}}\dfrac{\Gamma(v_1+\frac12)}{v_2^{v_1+\frac12}}
\left(1 -\dfrac{v_3^2}{2v_2v_4}\right)^{-(v_1+\frac12)}
\end{align*}
置いたとき, 事前分布に対する予測分布を
\begin{align*}
p(X_\ast|v)
&= \dfrac{1}{\sqrt{2\pi}}\sqrt{\dfrac{v_4}{v_4+1}}\dfrac{\Gamma(v_1+1)}{\Gamma(v_1+\frac12)}
\left(v_2 -\dfrac{v_3^2}{2v_4}\right)^{-\frac12}\\
&\quad\times \left\{ 1+ \dfrac{v_4}{2(v_4+1)(v_2-\frac{v_3^2}{2v_4})}\left(X_\ast-\frac{v_3}{v_4}\right)^2\right\}^{-(v_1+1)}
\end{align*}
と得ました.
今回記事では, ハイパーパラメータを変数変換することにより, この意味を明確にしていきましょう.
まず, 事前分布は次のようにガンマ分布の確率密度関数 と, 正規分布の確率密度関数 との積に書き換えられます.
\begin{align*}
p(\lambda|v )
&= z(v)\sqrt{\dfrac{2\pi}{v_4}} \lambda^{v_1-\frac12}e^{-(v_2-\frac{v_3^2}{2v_4})\lambda}
\exp\left\{-\dfrac{v_4\lambda}{2}\left(\mu - \frac{v_3}{v_4}\right)^2 \right\}\\
&= \textrm{G}\left(\lambda|v_1+\frac12, v_2-\frac{v_3^2}{2v_4}\right)
N\left(\mu| \frac{v_3}{v_4}, (v_4\lambda)^{-1}\right).
\end{align*}
このことから, 共役事前分布はハイパーパラメータを元にガンマ分布に従う分散を求め,
その分散も用いて正規分布に従う平均を求めるグラフィカルモデルと位置づけることが出来ます.
それでは, 事前分布が
\begin{equation}
p(\lambda|v )=\textrm{G}\left(\lambda|a, b\right)
N\left(\mu| m, (\beta\lambda)^{-1}\right)
\end{equation}
となるようなハイパーパラメータ()に置き換えましょう;
\begin{equation}
\left\{\begin{split}
&a=v_1+\frac12 \\
&b=v_2 - \frac{v_3^2}{2v_4} \\
&m=\frac{v_3}{v_4}\\
&(\beta\lambda)^{-1} = (v_4\lambda)^{-1}
\end{split}\right.
\Rightarrow
\left\{\begin{split}
&v_1= a-\frac12\\
&v_2=b +\frac12 m^2\beta\\
&v_3 = m\beta\\
&v_4= \beta
\end{split}\right.
\end{equation}
すると, 事前分布による予測分布は,
\begin{align*}
p(X_\ast|v)
&=\dfrac{\beta}{2\pi(\beta+1)b}\dfrac{\Gamma(a+\frac12)}{\Gamma(a)}\\
&\quad\times \left\{ 1+\dfrac{\beta}{2(\beta+1)b}(X_\ast-m)^2\right\}^{-(a+\frac12)}
\end{align*}
という, 自由度 の -分布の確率密度関数に近い形をしています.
実はこれは, 自由度 の -分布に従う確率変数 に対して,
\begin{equation}
\tilde{T}=m+\dfrac{\beta a}{(\beta+1)b}T
\end{equation}
により定められる確率変数 の確率密度関数です.
さらに, 事後分布 による予測分布 を求めましょう.
の代わりに を代入すればよいので,
新たにパラメータ を次のように置きます;
\begin{equation}
\left\{\begin{split}
&\hat{a}=(v+\hat{v})_1+\frac12=a+\dfrac{N}{2} \\
&\hat{b}=(v+\hat{v})_2 - \frac{(v+\hat{v})_3^2}{2(v+\hat{v})_4}
=b+\dfrac12\left(\sum X_n^2 +m^2\beta -\hat{m}^2\hat{\beta}\right) \\
&\hat{m}=\frac{(v+\hat{v})_3}{(v+\hat{v})_4}=\dfrac{m\beta+\sum X_n}{\beta+N}\\
&\hat{\beta} = (v+\hat{v})_4=\beta+N
\end{split}\right.
\end{equation}
すると, 事後分布による予測分布は,
\begin{align*}
p(X_\ast|\bar{X}, v)
&=\dfrac{\hat{\beta}}{2\pi(\hat{\beta}+1)\hat{b}}\dfrac{\Gamma(\hat{a}+\frac12)}{\Gamma(\hat{a})}\\
&\quad\times \left\{ 1+\dfrac{\hat{\beta}}{2(\hat{\beta}+1)\hat{b}}(X_\ast-\hat{m})^2\right\}^{-(\hat{a}+\frac12)}
\end{align*}
です.
したがって, 事後分布による予測分布はスチューデントの -分布に従い, 自由度, 中心, スケールパラメータはそれぞれ
\begin{align*}
2\hat{a} &= 2a + N\\
&= N + O(1), \quad \text{as } N\rightarrow \infty,\\
\hat{m} &= \dfrac{m\beta+\sum X_n}{\beta+N} \\
&= \dfrac{1}{N}\sum_{n=1}^N X_n +o(1), \quad \text{as } N\rightarrow \infty,\\
\dfrac{\hat{\beta}}{2(\hat{\beta}+1)\hat{b}}
&=\dfrac{N}{(N+\beta)(2b+\sum X_n^2 +m^2\beta -\hat{m}^2\beta)}\\
&= \dfrac{1}{N} \dfrac{N-1}{\displaystyle \sum_{n=1}^N\left( X_n -\frac{1}{N} \sum_{k=1}^N X_k\right)^2}+o(1), \quad \text{as } N\rightarrow \infty
\end{align*}
です.
, はそれぞれ平均, 分散の一致推定量なので, 事後分布は真の平均, 真の分散を持つような正規分布に漸近することが予想できます.

【ベイズ推論】一次元ガウス分布の平均・分散についてのベイズ推論の公式

正規分布についてのベイズ推論を考えましょう. 母数が つになることに気をつけましょう. 平均 , 精度 がともに未知の一次元ガウス分布を考えます.
\begin{equation}
N(x|\mu,\lambda^{-1})
= \sqrt{\dfrac{1}{2\pi}}\exp\left\{ \dfrac12\ln{\lambda} -\dfrac{\lambda}{2}x^2 + \lambda\mu x -\dfrac{1}{2}\lambda\mu^2\right\}
\end{equation}
により,
\begin{align*}
f(x)&=\sqrt{\dfrac{1}{2\pi}}, \\
\phi_1(\mu,\lambda)&= \ln{\lambda},\quad \phi_2(\mu,\lambda) =-\lambda,\quad
\phi_3(\mu,\lambda)= \lambda\mu,\quad \phi_4(\mu,\lambda) =-\frac12\lambda\mu^2,\\
g_1(x) &= \dfrac12, \quad g_2(x)=\dfrac12 x^2,\quad
g_3(x)=x,\quad g_4(x)=1
\end{align*}

なる指数型分布族です.

【ベイズ推論】指数型分布族を用いた共役事前分布の導出と予測分布の公式について

ベイズ推論をします. パラメータ をもつような確率変数 に対する独立同分布な族 を考えます. 事後分布を書く際にはハイパーパラメータ をしばしば省略します. すなわち とかきましょう.

指数型分布族と共役事前分布について

指数型分布族の定義

実数値確率変数 が未知パラメータ をもつ(条件付き)確率分布 に従うとします(つまり, 観測モデルを と定めます).
このとき, が指数型分布族に属するとは, ある値域の次元が等しい関数 , により,
\begin{align*}
p(x| \theta)
&= f(x)\exp\langle \phi(\theta) , g(x)\rangle ,\\
f(x)
&= \left(\int \exp\langle \phi(\theta) , g(x)\rangle\, d\theta\right)^{-1}
\end{align*}
が成立することをいいます. ガウス分布やポアソン分布などの代表的な分布はこの定義に属し, 具体例については次回また書いていきます. 今回は一般論として, この性質だけからベイズ推論における共役事前分布が得られ, さらにそれを用いた予測分布の公式が得られることを示していきます.

Python で1次元の熱方程式を2通りの方法で解いてみた

両端が断熱的な導線の熱伝導を記述する, 次の斉次 Neumann 境界付きの熱方程式の初期値境界値問題

を考えます*1.
ここで, , を表します.
前進差分法とスペクトル法を用いる方法の2通りでの数値計算例のデモを作成しました. *2
コードがメインで, 数値解析の手法については多くを説明しません.
意味が通らない箇所があれば気軽にコメントをよろしくお願いします.

*1: の第1引数を時間変数, 第2引数を空間変数としています.

*2:有限要素法でもやりたかった

不偏標本分散の一致性を示してみる

本記事では省略されがちであったり, 正規分布であることを仮定してあったりすることの多い, 不偏標本分散の母分散の推定量としての一致性を証明していきます.

Rにたくさんさいころを振らせて、さいころの分散を求めてみる

前回記事Rにさいころをたくさん振らせて、さいころの平均を求めてみる - わかばめにっきでは、Rの擬似乱数を使って標本平均が母平均の一致推定量であることを示しました。今回は後で述べる不偏標本分散という量が母分散の不偏推定量であることを調べていきます。