第一種変形ベッセル関数が方向統計学に登場する理由
ビルドアップ
変形ベッセル関数
$$ J_{\nu}(x) = \sum \limits_{n=0}^{\infty} \frac{(-1)^{n} }{\Gamma (n+1) \Gamma (n+\nu+1)} \left(\frac{x}{2} \right)^{2n+\nu} $$
第一種ベッセル関数 $J_{\nu}$ に対して次のように定義される $I_{\nu}$ を変形第一種ベッセル関数という1。
$$ \begin{align*} I_{\nu} (z) :=& i^{-\nu} J_{\nu} \left( iz \right) \\ =& \left( {{ z } \over { 2 }} \right)^{\nu} \sum_{k=0}^{\infty} {{ {{ z } \over { 2 }}^{2k} } \over { k! \Gamma \left( \nu + k + 1 \right) }} \\ =& {{ \left( {{ z } \over { 2 }} \right)^{\nu} } \over { \sqrt{\pi} \Gamma \left( \nu + {{ 1 } \over { 2 }} \right) }} \int_{-1}^{1} e^{zt} \left( 1 - t^{2} \right)^{\nu - {{ 1 } \over { 2 }}} dt \end{align*} $$
方向統計
一方、方向統計学statisticalとは、通常のユークリッド空間ではなく、ある多様体上の確率分布や統計的推測などを研究する分野だ。たとえば地球のような球体を代表するスフィアや、$2 \pi$-モジュラスで表されるトーラスなどに置かれたデータを扱うが、すぐに空間情報統計学(球面)や分⾚間の角度(トーラス)に応用できるなど、その未来は明るい部門だ。だが、ここに登場する分布は、次のような奇妙な確率密度関数を持っている。 $$ f_{p} \left( \mathbf{x} ; \mu , \kappa \right) := \left( {{ \kappa } \over { 2 }} \right)^{p/2 - 1} {{ 1 } \over { \Gamma \left( p/2 \right) I_{p/2 - 1} \left( \kappa \right) }} \exp \left( \kappa \mu^{T} \mathbf{x} \right) \qquad , \mathbf{x} \in S^{p-1} $$ ここで最初の複雑な因子は、修正第一種ベッセル関数 $I_{\nu}$ を含むが、$\int_{S^{p-1}} f d \mathbf{x} = 1$ になるように定数で正規化normalizeしてくれるものだ。このように複雑な関数が使用される理由は単純だ。
解答
第一種ベッセル関数の導出: $\nu \in \mathbb{R}$ に対して、下のような形の微分方程式を$\nu$次ベッセル方程式 という。 $$ \begin{align*} && x^{2} y^{\prime \prime} +xy^{\prime}+(x^{2}-\nu^{2})y &= 0 \\ \text{or} && y^{\prime \prime}+\frac{1}{x} y^{\prime} + \left( 1-\frac{\nu^{2}}{x^{2}} \right)y &= 0 \end{align*} $$ ベッセル方程式は、球面座標系で波動方程式を解く時に登場する微分方程式だ。 係数は定数ではなく、独立変数 $x$ に依存する。フロベニウスの方法で解を求めることができ、級数解 は次のようになる。 $$ \begin{align*} J_{\nu}(x) &= \sum \limits_{n=0}^{\infty} \frac{(-1)^{n} }{\Gamma (n+1) \Gamma (n+\nu+1)} \left(\frac{x}{2} \right)^{2n+\nu} \\ J_{-\nu}(x) & =\sum \limits_{n=0}^{\infty}\frac{(-1)^{n}}{\Gamma (n+1)\Gamma (n-\nu+1)} \left( \frac{x}{2} \right)^{2n-\nu} \end{align*} $$
ベッセル関数の導出を考えると、ベッセル方程式 という微分方程式の一つやその解など、難しい表現が登場する。しかし、方向統計学に関する上記の引用で重要な文はただ一つだ。
“ベッセル方程式は、球面座標系で波動方程式を解く時に登場する微分方程式だ。”
数理物理学で波動方程式について語るかもしれないが、我々に必要なのは$\int_{S^{p-1}} f d \mathbf{x} = 1$だけだ。問題は、通常のユークリッド空間と違って球面上の確率密度関数の値が「中心から離れて$0$に近づく」という現象がないため、積分自体が簡単ではないことだ。正規分布の確率密度関数が円$S^{1}$を包む形を想像してみてほしい。
上の図で$\tau = 1$の時を見ると、正規分布の果てしない尻尾は、無限に$S^{1}$を回りながら、無限に$0$ではない厚さを加えている。2これは円周率の二倍である$2 \pi$を周期として繰り返し、これがまさに「球面上の波」と同じ形を持つ理由であり、ベッセル関数が使用できる理由だ。
Sungkyu Jung. “Geodesic projection of the von Mises–Fisher distribution for projection pursuit of directional data.” Electron. J. Statist. 15 (1) 984 - 1033, 2021. https://doi.org/10.1214/21-EJS1807 ↩︎
Straub, J. (2017). Bayesian Inference with the von-Mises-Fisher Distribution in 3D. https://www.semanticscholar.org/paper/Bayesian-Inference-with-the-von-Mises-Fisher-in-3-D-Straub/26d5bb31153df418388b6eb242b2d8842c039c2d#extracted ↩︎