Von Mises分布の確率密度関数を導出してみた。

前置き

この記事で整数次の(第1種)変形Bessel関数からなる関数列の母関数をTaylor展開すると、Von Mises分布の展開式を導出できて、かつそれを使うとVon Mises分布の累積分布関数を求めることができることを示しました。

また、「Von Mises分布の確率密度関数は2変量の正規分布の確率密度関数から導出できる。」ともののWikipedia等に書いてありますが、具体的な方法については書いていなかったりします。

そこで、具体的に計算してみることにしました。

計算してみます。

2変量正規分布の確率密度関数

ここで唐突ですが、平均$\boldsymbol{\mu}$、正定値の分散共分散行列$\Sigma = \sigma^2 I = \left( \begin{array}{cc} \sigma^2 & 0 \cr 0 & \sigma^2 \end{array}\right)$($I$は$2\times 2$の単位行列)、確率変数を$\boldsymbol{x}$とすると、2変量正規分布の確率密度関数$p(\boldsymbol{x})$は… \begin{align} p(\boldsymbol{x}) &= \frac{1}{2\pi\sqrt{\det\Sigma}}\exp\left[ -\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T\Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu}) \right] \nonumber\cr &= \frac{1}{2\pi\sigma^2}\exp\left[ -\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T\left( \begin{array}{cc} \displaystyle\frac{1}{\sigma^2} & 0 \cr 0 & \dfrac{1}{\sigma^2} \end{array} \right)(\boldsymbol{x}-\boldsymbol{\mu}) \right] \nonumber\cr &= \frac{1}{2\pi\sigma^2}\exp\left[ -\frac{1}{2\sigma^2}(\boldsymbol{x}-\boldsymbol{\mu})^2\right] \nonumber\cr &= \frac{1}{2\pi\sigma^2}\exp\left[ -\frac{1}{2\sigma^2}(\|\boldsymbol{x}\|^2+\|\boldsymbol{\mu}\|^2-2\|\boldsymbol{x}\|\|\boldsymbol{\mu}\|\cos\theta)\right] \label{eq:mnddef} \end{align} と表すことができます。ただし、$\theta$は$\|\boldsymbol{x}\|$と$\|\boldsymbol{\mu}\|$が2次元の直交座標系(以下、単に「座標平面」と書きます。)においてなす角を表します。

ここまでに登場した数で座標平面に書くことができるもので、$\mu(=\phi-\theta)$以外のものを座標平面上に図示すると以下のような感じになります。

変数変換

確率変数$\boldsymbol{x}$が持つ情報のうち、Von Mises分布では座標平面上$x軸$正の向きと$\boldsymbol{x}$のなす角($\phi$とします。)のみが必要ですので、(\ref{eq:mnddef})式を$\phi$を使った形に変形できないか考えます。そこで、$\boldsymbol{x} = (r\cos\phi,r\sin\phi)^T$とおくと…

\begin{align} p(\boldsymbol{x})d\boldsymbol{x} &= p(\boldsymbol{x}(r,\phi))\left|\frac{\partial\boldsymbol{x}}{\partial r\partial\phi}\right|drd\phi \nonumber\cr &= p(\boldsymbol{x}(r,\phi))\left| \begin{array}{cc} \cos\phi & -r\sin\phi \cr \sin\phi & r\cos\phi \end{array} \right|drd\phi \nonumber\cr &= p(\boldsymbol{x}(r,\phi))rdrd\phi \label{eq:probconv} \end{align}

と変換できます。さらに、座標平面上$x軸$正の向きと$\boldsymbol{\mu}$のなす角を$\mu$とおくと、$\theta=\phi - \mu$となることに留意しつつ、(\ref{eq:probconv})式に(\ref{eq:mnddef})式を適用すると…

\begin{align} p(\boldsymbol{x}(r,\phi))r &= \frac{1}{2\pi\sigma^2}r\exp\left[ -\frac{1}{2\sigma^2}(r^2+\|\boldsymbol{\mu}\|^2-2r\|\boldsymbol{\mu}\|\cos(\phi - \mu))\right] \label{eq:mnddefsecond} \end{align}

となります。

定数とみなす部分の置き換え。

ここからは、角度($\phi$)以外の数は適当に置き換えつつ式変形します。

まず、(\ref{eq:mnddefsecond})式を

\begin{align} p(\boldsymbol{x}(r,\phi))r &= \frac{r}{2\pi\sigma^2}\exp\left[ -\frac{1}{2\sigma^2}(r^2+\|\boldsymbol{\mu}\|^2)\right]\exp\left[-\frac{1}{2\sigma^2}(-2r\|\boldsymbol{\mu}\|\cos(\phi - \mu))\right] \nonumber\cr &= \frac{r}{2\pi\sigma^2}\exp\left[ -\frac{1}{2\sigma^2}(r^2+\|\boldsymbol{\mu}\|^2)\right]\exp\left[\frac{r\|\boldsymbol{\mu}\|\cos(\phi - \mu)}{\sigma^2}\right]\label{eq:mnddefthird} \end{align}

のように変形し、$\dfrac{r}{2\pi\sigma^2}\exp\left[ -\dfrac{1}{2\sigma^2}(r^2+\|\boldsymbol{\mu}\|^2)\right] = C$とおきます。また、(\ref{eq:mnddefthird})式の左辺についてもあらためて$p(\phi)$とおき直します。

すると、(\ref{eq:mnddefthird})式は…

\begin{align} p(\phi) &= C\exp\left[\frac{r\|\boldsymbol{\mu}\|\cos(\phi - \mu)}{\sigma^2}\right]\label{eq:mnddeffourth} \end{align}

と書き換えることができます。

次に$\exp$の指数部の$\phi$とは無関係な部分(=定数とみなせる)を

\begin{align} \frac{r\|\boldsymbol{\mu}\|}{\sigma^2} &= m \label{eq:m} \end{align}

とおくと、(\ref{eq:mnddeffifth})式のように変形できます。

\begin{align} p(\phi) &= C\exp\left[m\cos(\phi - \mu)\right]\label{eq:mnddeffifth} \end{align}

確率密度関数の導出。

(\ref{eq:mnddeffifth})式より、$p(\phi)$は$0 \le \phi \le 2\pi$において負でない値をとり、かつ$p(\phi+2\pi) = p(\phi)$になることから、

\begin{align} \int_0^{2\pi} p(\phi)d\phi = 1 \label{eq:density} \end{align}

となるように$C$を定めると、$p(\phi)$を$\phi$を確率変数とする確率密度関数とすることができます。

そこで、(\ref{eq:density})式に(\ref{eq:mnddeffifth})式を代入して計算します。

すると…

\begin{align} \int_0^{2\pi} p(\phi)d\phi &= \int_0^{2\pi} C\exp\left[m\cos(\phi - \mu)\right] d\phi \nonumber\cr &= \int_{-\mu}^{2\pi-\mu} C\exp\left[m\cos t\right] dt \nonumber\cr &= \int_{0}^{2\pi} C\exp\left[m\cos t\right] dt\nonumber\cr &= \int_{0}^{\pi} C\exp\left[m\cos t\right] dt + \int_{\pi}^{2\pi} C\exp\left[m\cos t\right] dt \label{eq:densitysecond} \end{align}

となります($p(\phi+2\pi) = p(\phi)$であることを用いています)。(\ref{eq:densitysecond})式の$t$を$\phi$と書き換えて、右辺第2項について$t = 2\pi - \phi$と置き換えると、$\cos(2\pi - \phi) = \cos\phi$となることから…

\begin{align} \int_0^{2\pi} p(\phi)d\phi &= \int_{0}^{\pi} C\exp\left[m\cos \phi\right] d\phi - \int_{\pi}^{0} C\exp\left[m\cos(2\pi - t)\right] dt \nonumber\cr &= \int_{0}^{\pi} C\exp\left[m\cos \phi\right] d\phi + \int_{0}^{\pi} C\exp\left[m\cos t\right] dt \nonumber\cr &= 2\int_{0}^{\pi} C\exp\left[m\cos \phi\right] d\phi\label{eq:densitythird} \end{align}

(\ref{eq:densitythird})式の左辺は初等的に解くことはできませんが、0次の第1種変形Bessel関数$I_0(m)$が

\begin{align} I_0(x) &= \frac{1}{\pi}\int_0^{\pi}exp(x\cos \theta)d\theta \label{eq:mbffzero} \end{align}

と書けることを利用すると、(\ref{eq:densitythird})式は…

\begin{align} \int_0^{2\pi} p(\phi)d\phi &= 2\pi CI_0(m) \label{eq:densityfourth} \end{align}

と書くことができます。これが1に等しくなるときの$C$の値を求めると、

\begin{align} C = \frac{1}{2\pi I_0(m)} \label{eq:densityC} \end{align}

となりますので、

\begin{align} p(\phi) &= \frac{\exp\left[m\cos(\phi - \mu)\right]}{2\pi I_0(m)}\label{eq:densityfinal} \end{align}

となります。

(\ref{eq:densityfinal})式がVon Mises分布の確率密度関数になります。$\blacksquare$

まとめ

Wikipediaの記事やPRMLなどを読んでいてもいまいち$x_1 = r\cos\theta, x_2 = r\sin\theta$とおいて確率変数を変換した後の$r$の処理の方法がいまいち理解できていなかったのですが、「確率変数の変数変換を行った後は、$r$を変数とみなさず固定して考える。」こととしました。

これで、Von Mises分布に従う乱数を生成するコード(以下、単に「コード」と書きます。)をscratchから実装できるような気がします。👍

(\ref{eq:m})式の$m$は集中度と呼ばれるパラメータですが、乱数生成時にこれと$\|\boldsymbol{\mu}\|$のみが与えられた場合には、$\sigma^2$と$r$の比を定める必要があります。おそらく、$r=1$とおくことにより$\sigma$の値を定めて、正規分布$N_2(\boldsymbol{\mu},\Sigma)$に従う乱数を2個発生させてその向きを求めれば良いのではないか… と思います(試していないのでよくわかりません。気が向いたら試してみます)。

この記事は以上です。

リンク

panda大学習帳外伝

素材の味を活かしたリール動画。