風向計が今日の朝、350度（ほぼ北）を指していました。昼には10度（やはりほぼ北）を指しました。この2つの方向の「平均」は何度でしょうか。算術平均で計算すると $(350 + 10)/2 = 180$ 度、つまり真南という答えが返ってきます。北寄りの2つの観測から平均が南になる——明らかに何かがおかしいわけです。
左パネルが示す通り、直線上の算術平均では 350° と 10° の中間である 180°（南）が答えとなり、物理的に全く意味をなしません。右パネルの円周上では、2つの観測を単位ベクトルに変換して合成すると、その方向は 0°（北）を正確に指します。これが「円周上の平均」の本質的な考え方です。
おかしいのは平均の取り方ではなく、角度に通常の正規分布を当てはめるという発想そのものです。角度は直線上の数ではなく、円周上を一周してつながる量で、$350° = -10°$ という同一視が成り立つ世界です。直線上の確率分布をそのまま持ち込めば、こうした0°/360°の継ぎ目で必ず破綻します。これを正面から解決するのが、円周 $S^1$ 上で定義された確率分布の代表選手、フォン・ミーゼス分布（von Mises distribution） です。
この分布は身近なところで広く活躍しています。たとえば、ロボットや自動運転車の方位角の推定（センサ観測値を確率的に統合する）、ヘディング推定、神経科学でニューロン発火の位相同期解析、地質学で測定された古地磁気の方向データの解析、画像処理ではエッジ方位のヒストグラム解析、深層学習では「角度の損失関数」として用いられる Cosine Similarity の理論的裏付け、いずれもフォン・ミーゼス分布が背後にいます。
本記事の内容
- 角度データに正規分布を当てはめると何が壊れるかの具体例
- 円周 $S^1$ 上の最大エントロピー分布としてフォン・ミーゼス分布を導出
- 第1種変形ベッセル関数 $I_0$ が正規化定数として現れる理由
- 集中度パラメータ $\kappa$ の意味、$\kappa \to 0$ で一様分布、$\kappa \to \infty$ で正規近似
- 最尤推定（$\mu$ は閉形式、$\kappa$ はベッセル関数比のニュートン法）
- 球面 $S^{n-1}$ 上のフォン・ミーゼス・フィッシャー分布への拡張
- Python での標本生成・極座標可視化・最尤推定・ベイズ姿勢推定への応用
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 正規分布とは？定義やグラフの見方を徹底解説
- 多変量正規分布の理論と実装
- ガウス分布の線形変換と再生性
- ベルヌーイ分布と最尤推定
- 剛体の姿勢表現（オイラー角・回転行列・クォータニオン）
- 姿勢決定アルゴリズム
直感 — なぜ角度に正規分布を当てはめてはいけないのか
冒頭の風向計の話を、もう少し丁寧に追ってみましょう。観測値 $\theta_1 = 350°$、$\theta_2 = 10°$ の2つを「真の方向 $\mu$ のまわりに正規分布で揺らぐ観測ノイズ」だと仮定すると、最尤推定は
$$ \hat{\mu} = \arg\max_\mu \prod_{i} \frac{1}{\sqrt{2\pi}\sigma}\exp\!\Big(-\frac{(\theta_i – \mu)^2}{2\sigma^2}\Big) = \frac{1}{2}(\theta_1 + \theta_2) = 180° $$
となります。最後の等号は正規分布の平均推定がサンプル平均になるという有名な事実です。しかし、$\theta_1 = 350°$ を $\theta_1 = -10°$ と書き換えて同じ計算をすれば $\hat{\mu} = 0°$ になります。同じ物理現象なのに、角度の表記を変えただけで答えが180度ずれる——これは推定が観測の表記に依存しているという致命的な欠陥です。
問題の根源は、正規分布の密度関数が距離 $(\theta_i – \mu)^2$ をユークリッド距離で測っていることです。円周上では、$10°$ と $350°$ の「真の」距離は20度（短い弧の方）であって340度ではありません。直線上の確率分布は、円周上の自然な距離 $|\theta_i – \mu|_{\bmod 2\pi} = \min(|\theta_i – \mu|,\, 2\pi – |\theta_i – \mu|)$ を表現できないのです。
ではどうするか。一案として「観測ごとに最も近い表記に正規化する」というアドホックな処置もありますが、これは確率モデルとしての一貫性を失います（積分が1にならない、ベイズ更新で爆発する、など）。本筋は 円周 $S^1$ 上で最初から定義された確率分布 を作ることです。次節では、$S^1$ 上で「最も自然な」正規分布の対応物を導出します。
4つのパネルは $\kappa = 0.5, 2, 5, 20$ のフォン・ミーゼス密度関数を示しています。左端の $\kappa=0.5$ はほぼ水平（一様分布に近い）で、右端の $\kappa=20$ では鋭い山が立ち、正規分布のベル型に見えます。これは $\kappa \to \infty$ で正規分布に漸近するという性質の視覚的確認です。緑の点線（正規近似）が $\kappa=5, 20$ で密度関数によく重なっていることも確認できます。
フォン・ミーゼス分布の定義と導出
円周上で「正規分布の代わり」になる分布
直線上で正規分布 $\mathcal{N}(\mu, \sigma^2)$ が特権的な地位を持つ理由はいくつかありますが、最も深いのは「平均と分散を固定したときに最大エントロピーを持つ分布である」という性質です。最大エントロピーとは「与えられた制約以外には何も決めつけない、最も無情報な分布」を意味します。
円周上で同じ発想を働かせます。角度 $\theta$ の代わりに、その単位ベクトル表現 $\bm{u} = (\cos\theta, \sin\theta)^\top$ を使えば、円周上の点は2次元平面上の単位円周上の点として扱えます。すると「方向 $\mu$ のまわりに集中している」という性質は、ベクトルの平均 $\langle \bm{u} \rangle$ が方向 $\bm{\mu} = (\cos\mu, \sin\mu)^\top$ にどれだけ揃っているか、つまり内積
$$ \langle \bm{u} \cdot \bm{\mu} \rangle = \langle \cos(\theta – \mu) \rangle $$
で測れます。この量を固定した上で、エントロピーを最大化する円周上の確率密度を求めるのです。
最大エントロピー導出
円周上の確率密度 $f(\theta)$ について、エントロピー
$$ H[f] = -\int_0^{2\pi} f(\theta)\log f(\theta)\,d\theta $$
を、次の制約のもとで最大化します。
$$ \int_0^{2\pi} f(\theta)\,d\theta = 1, \qquad \int_0^{2\pi} \cos(\theta – \mu)\,f(\theta)\,d\theta = R $$
第1式は確率の正規化、第2式は「平均的に方向 $\mu$ を向いている度合い」を $R$ に固定するものです。ラグランジュ未定乗数法で
$$ \mathcal{L}[f] = -\int f\log f\,d\theta – \lambda_0\!\left(\int f\,d\theta – 1\right) – \lambda_1\!\left(\int \cos(\theta – \mu)f\,d\theta – R\right) $$
を作り、$f$ について変分を取ります。被積分関数を $f$ で微分するイメージで $\delta\mathcal{L}/\delta f = 0$ を解くと
$$ -\log f(\theta) – 1 – \lambda_0 – \lambda_1 \cos(\theta – \mu) = 0 $$
両辺の指数を取ると、正規化定数を一つにまとめて
$$ f(\theta) = C\,\exp\!\big(\kappa \cos(\theta – \mu)\big) $$
の形になります（$\kappa = -\lambda_1$ と置き換えました）。これがフォン・ミーゼス分布の密度関数です。最大エントロピー的に「平均方向 $\mu$ と集中度 $R$ 以外には何も決めつけない」という条件から、$\cos(\theta – \mu)$ の指数関数という形が一意に導かれることがわかります。
正式な定義
フォン・ミーゼス分布の密度関数を、正規化定数を明示して書くと
$$ \boxed{\,f(\theta\mid\mu,\kappa) = \frac{1}{2\pi I_0(\kappa)}\exp\!\big(\kappa\cos(\theta – \mu)\big), \quad \theta \in [0, 2\pi)\,} $$
となります。$\mu \in [0, 2\pi)$ は平均方向（モード位置）、$\kappa \geq 0$ は集中度パラメータ（concentration、正規分布における $1/\sigma^2$ のような役割）、$I_0(\kappa)$ は次節で詳しく扱う第1種変形ベッセル関数 で、円周上で積分を 1 に正規化する役目を持ちます。
定義式の構造を読み解くと、$\cos(\theta – \mu)$ は $\theta = \mu$ で最大値1、$\theta = \mu \pm \pi$（反対側）で最小値 $-1$ を取ります。指数関数を介すことで、$\theta = \mu$ 付近に確率質量が集中し、反対側で最小になります。$\kappa$ が大きいほど指数関数の感度が増し、ピークがシャープになります。次節では、この見かけの構造を保証する正規化定数 $I_0(\kappa)$ がなぜこの形になるのかを丁寧に追います。
方向データは極座標で可視化するのが最も自然です。左パネル（$\kappa=4$）では北東（45°）方向に強い集中が見られ、右パネル（$\kappa=1.5$）では南西（200°）寄りながらもかなり広がりがあることが一目でわかります。直交座標で同じデータを描くと 0°/360° の継ぎ目で分断されますが、極座標なら方向の連続性がそのまま現れます。
ベッセル関数と正規化定数
$I_0(\kappa)$ の定義
第1種変形ベッセル関数（modified Bessel function of the first kind）$I_n(\kappa)$ は、もともと変形ベッセル微分方程式 $x^2 y” + x y’ – (x^2 + n^2)y = 0$ の解として登場しますが、ここでは積分表示
$$ I_n(\kappa) = \frac{1}{\pi}\int_0^\pi e^{\kappa\cos\phi}\cos(n\phi)\,d\phi $$
が便利です。$n = 0$ の場合に注目すると
$$ I_0(\kappa) = \frac{1}{\pi}\int_0^\pi e^{\kappa\cos\phi}\,d\phi = \frac{1}{2\pi}\int_0^{2\pi} e^{\kappa\cos\phi}\,d\phi $$
となります（被積分関数が $\phi$ について偶対称なので積分区間を倍にできます）。一番右の式は、まさに我々のフォン・ミーゼス密度の正規化に必要な積分の形そのものです。
正規化定数の導出
$f(\theta\mid\mu,\kappa) = C\exp(\kappa\cos(\theta – \mu))$ を $[0, 2\pi)$ で積分し、1になるように $C$ を決めます。変数変換 $\phi = \theta – \mu$ を施すと（周期関数の積分なので積分区間はそのまま $[0, 2\pi)$ で済みます）
$$ \int_0^{2\pi} C\exp(\kappa\cos(\theta – \mu))\,d\theta = C\int_0^{2\pi} e^{\kappa\cos\phi}\,d\phi $$
最右辺の積分は前項の $I_0(\kappa)$ の積分表示そのものなので、$\int_0^{2\pi} e^{\kappa\cos\phi}d\phi = 2\pi I_0(\kappa)$ が直ちに得られます。したがって正規化条件 $C \cdot 2\pi I_0(\kappa) = 1$ から
$$ C = \frac{1}{2\pi I_0(\kappa)} $$
となり、冒頭で書いたフォン・ミーゼス分布の密度の係数が確定します。$I_0$ という見慣れない関数が出てきたのは、$e^{\kappa\cos\phi}$ という積分が初等関数では閉形式で書けないからで、ベッセル関数はその「名前」を与えただけ、と理解すると気が楽です。
$I_0(\kappa)$ の振る舞い
数値計算では scipy.special.i0(kappa) で簡単に評価できますが、極限の振る舞いだけ押さえておくと感覚がつかめます。級数展開は
$$ I_0(\kappa) = \sum_{k=0}^{\infty}\frac{(\kappa/2)^{2k}}{(k!)^2} = 1 + \frac{\kappa^2}{4} + \frac{\kappa^4}{64} + \cdots $$
で、$\kappa \to 0$ で $I_0(\kappa) \to 1$。一方、$\kappa \to \infty$ では漸近展開
$$ I_0(\kappa) \sim \frac{e^\kappa}{\sqrt{2\pi\kappa}}\left(1 + \frac{1}{8\kappa} + \cdots\right) $$
で発散します。後で出てくる近似評価でこの2つの極限を使います。次節では、この $I_0(\kappa)$ を通して集中度パラメータ $\kappa$ の意味をさらに掘り下げ、$\kappa \to 0$ と $\kappa \to \infty$ の極限でフォン・ミーゼス分布が何に化けるかを見ます。
左パネルでは $I_0(\kappa)$ と $I_1(\kappa)$ が単調増加し、大きな $\kappa$ では漸近式 $e^\kappa/\sqrt{2\pi\kappa}$（黄点線）に接近していくことがわかります。右パネルの片対数表示では、指数関数的な成長が直線として現れ、$I_0$ が $I_1$ より常に大きい（$I_0 > I_1 > 0$）ことが読み取れます。この大小関係は $A(\kappa) = I_1/I_0 < 1$ という性質、つまり平均合成ベクトル長が必ず 1 未満になることを保証します。
集中度パラメータ $\kappa$ と二つの極限
$\kappa \to 0$ — 一様分布への退化
$\kappa = 0$ の場合、密度関数は
$$ f(\theta\mid\mu, 0) = \frac{1}{2\pi I_0(0)}\exp(0\cdot\cos(\theta – \mu)) = \frac{1}{2\pi \cdot 1}\cdot 1 = \frac{1}{2\pi} $$
となり、$\theta$ によらず一定です。つまり円周上の一様分布に一致します。これは「方向に関する情報がない」状態に対応します。観測がまったく無情報になる場合や、ベイズ姿勢推定で事前分布として「方向についてはノーアイデア」と表したい場合に、$\kappa = 0$ のフォン・ミーゼス分布が自然な選択肢になります。
3パネルの比較で両極限が一目瞭然です。$\kappa=0.01$（左）ではフォン・ミーゼス密度が $1/(2\pi)$ の水平線（一様分布）と重なり、$\kappa=30$（右）では正規分布（緑点線）と見分けがつかないほど一致します。$\kappa=3$（中央）は両極端の中間で、「方向の好みはあるが、かなりばらつきもある」状態です。この視覚的比較によって $\kappa$ が「情報量の強さ」を制御するパラメータであることが直感的につかめます。
$\kappa \to \infty$ — 正規分布への漸近
逆に $\kappa$ が大きい極限を考えます。$\cos(\theta – \mu)$ をテイラー展開すると、$\theta \approx \mu$ の近傍で
$$ \cos(\theta – \mu) = 1 – \frac{(\theta – \mu)^2}{2} + \mathcal{O}((\theta – \mu)^4) $$
なので、密度は
$$ f(\theta\mid\mu,\kappa) \approx \frac{e^\kappa}{2\pi I_0(\kappa)}\exp\!\left(-\frac{\kappa(\theta – \mu)^2}{2}\right) $$
と書けます。係数 $e^\kappa /(2\pi I_0(\kappa))$ を $I_0$ の漸近展開 $I_0(\kappa) \sim e^\kappa/\sqrt{2\pi\kappa}$ で評価すると
$$ \frac{e^\kappa}{2\pi I_0(\kappa)} \approx \frac{e^\kappa}{2\pi \cdot e^\kappa/\sqrt{2\pi\kappa}} = \frac{\sqrt{2\pi\kappa}}{2\pi} = \frac{1}{\sqrt{2\pi/\kappa}} $$
となり、結果として
$$ f(\theta\mid\mu,\kappa) \approx \frac{1}{\sqrt{2\pi\sigma^2}}\exp\!\left(-\frac{(\theta-\mu)^2}{2\sigma^2}\right), \quad \sigma^2 = \frac{1}{\kappa} $$
という平均 $\mu$、分散 $1/\kappa$ の正規分布に漸近することがわかります。$\kappa$ が大きいということは「ピークが鋭く、$2\pi$ 周期の繰り返しが効かないほど局所化している」状態で、その極限では円周であることを忘れて直線上の正規分布として扱えるわけです。実用上は $\kappa \gtrsim 10$ あたりでこの近似が良く効きます。
平均合成ベクトルと $A(\kappa)$
集中度を直感的に量る量として、平均合成ベクトル長
$$ \rho = \mathbb{E}[\cos(\theta – \mu)] = \frac{1}{2\pi I_0(\kappa)}\int_0^{2\pi}\cos(\theta – \mu)\,e^{\kappa\cos(\theta-\mu)}\,d\theta $$
を定義します。$\partial/\partial\kappa$ を被積分関数の中に入れて計算すると（積分と微分の交換は可）
$$ \rho = \frac{1}{I_0(\kappa)}\cdot \frac{d I_0(\kappa)}{d\kappa} = \frac{I_1(\kappa)}{I_0(\kappa)} \equiv A(\kappa) $$
となります。ここで $I_1$ は $n=1$ の変形ベッセル関数。$A(\kappa) = I_1(\kappa)/I_0(\kappa)$ は$0$ から $1$ の間の単調増加関数で、$\kappa = 0$ で $A=0$（一様）、$\kappa \to \infty$ で $A \to 1$（デルタ関数化）になります。$A(\kappa)$ は次節の最尤推定で中心的役割を果たします。次節では、観測データから $\mu$ と $\kappa$ を推定する具体的な方法に踏み込みます。
左パネルは $A(\kappa) = I_1(\kappa)/I_0(\kappa)$ が 0 から 1 の間を単調に増加することを示しています。$\kappa=0.5$ では $A \approx 0.24$（まだかなり広がりがある）、$\kappa=5$ では $A \approx 0.79$（かなり集中）と読み取れます。右パネルは逆関数の視点で、観測から得た平均合成ベクトル長 $\bar{R}$ が 0.9 を超えると $\hat{\kappa}$ が急激に大きくなることを示しており、この非線形性がニュートン法による数値解法を必要とする理由です。
最尤推定と球面拡張（vMF）
対数尤度
$N$ 個の独立な観測 $\theta_1, \dots, \theta_N$ がフォン・ミーゼス分布から得られたと仮定します。対数尤度は
$$ \mathcal{L}(\mu, \kappa) = \sum_{i=1}^N \log f(\theta_i\mid\mu,\kappa) = -N\log(2\pi I_0(\kappa)) + \kappa\sum_{i=1}^N \cos(\theta_i – \mu) $$
となります。$\mu$ と $\kappa$ について順に最大化します。
$\mu$ の最尤推定 — 平均合成ベクトル
$\mu$ について偏微分してゼロと置きます。
$$ \frac{\partial \mathcal{L}}{\partial \mu} = \kappa\sum_{i=1}^N \sin(\theta_i – \mu) = 0 $$
$\sin(\theta_i – \mu) = \sin\theta_i\cos\mu – \cos\theta_i\sin\mu$ と展開して整理すると
$$ \cos\mu\sum_i \sin\theta_i – \sin\mu\sum_i\cos\theta_i = 0 $$
これは
$$ \tan\hat\mu = \frac{\sum_i \sin\theta_i}{\sum_i \cos\theta_i} $$
を意味します。象限の取り違えを防ぐため、実装では np.arctan2 を使うのが定番です。
$$ \boxed{\;\hat\mu = \mathrm{atan2}\!\left(\sum_i \sin\theta_i,\;\sum_i \cos\theta_i\right)\;} $$
この推定量の意味を視覚化すると、各観測 $\theta_i$ を単位ベクトル $\bm{u}_i = (\cos\theta_i, \sin\theta_i)^\top$ に変換してベクトル和を取り、その方向角を取る操作になります。冒頭の風向例 $\theta_1=350°, \theta_2=10°$ で計算すると、$\bm{u}_1 + \bm{u}_2 = (\cos350° + \cos10°,\;\sin350°+\sin10°) \approx (1.97, 0)$ から $\hat\mu = 0°$ が得られ、直感的に正しい「北」が出ます。直線上の算術平均が破綻していた問題が、ここでは自然に解消されています。
$\kappa$ の最尤推定 — ベッセル関数比の方程式
$\kappa$ についての最尤方程式は
$$ \frac{\partial \mathcal{L}}{\partial \kappa} = -N\frac{I_0′(\kappa)}{I_0(\kappa)} + \sum_i \cos(\theta_i – \hat\mu) = 0 $$
$I_0′(\kappa) = I_1(\kappa)$ なので、$\bar{R} \equiv \frac{1}{N}\sum_i \cos(\theta_i – \hat\mu)$ と書けば
$$ \boxed{\;A(\hat\kappa) \equiv \frac{I_1(\hat\kappa)}{I_0(\hat\kappa)} = \bar{R}\;} $$
という関係になります。$\bar{R}$ は 平均合成ベクトル長で、$\hat\mu$ 方向への観測の集中度を直接表す統計量です。$A(\kappa)$ が解析的に逆関数を持たないので、$\hat\kappa$ は数値的に求める必要があります。
Newton 法による $\hat\kappa$ の求解
$g(\kappa) \equiv A(\kappa) – \bar R = 0$ を解くニュートン法は
$$ \kappa_{k+1} = \kappa_k – \frac{A(\kappa_k) – \bar R}{A'(\kappa_k)} $$
です。導関数 $A'(\kappa)$ は、ベッセル関数の漸化式 $I_0′ = I_1$, $I_1′ = I_0 – I_1/\kappa$ から
$$ A'(\kappa) = \frac{I_1’I_0 – I_1 I_0′}{I_0^2} = \frac{(I_0 – I_1/\kappa)I_0 – I_1^2}{I_0^2} = 1 – \frac{A(\kappa)}{\kappa} – A(\kappa)^2 $$
と閉形式で得られます。初期値には、よく使われる近似式
$$ \hat\kappa_0 \approx \begin{cases} \bar R(2 – \bar R^2)/(1 – \bar R^2), & \bar R < 0.85 \\ 1/\big(2(1 - \bar R) - (1 - \bar R)^2 - (1 - \bar R)^3\big), & \bar R \geq 0.85 \end{cases} $$
を使うと、数回の反復で収束します。$\bar R$ が 1 に近い（観測がきれいに揃っている）ほど $\hat\kappa$ は大きくなり、$\bar R \to 0$ では $\hat\kappa \to 0$（観測がバラバラ → 集中度ゼロ → 一様分布）となります。
球面への拡張 — フォン・ミーゼス・フィッシャー分布
3次元以上の方向データを扱いたい場合（たとえばロボットの姿勢軸を3次元単位ベクトルで表す、テキスト埋め込みベクトルの方向、地磁気の3次元単位ベクトル）、円周 $S^1$ ではなく $n-1$ 次元球面 $S^{n-1}$ 上の分布が必要です。これがフォン・ミーゼス・フィッシャー分布（von Mises-Fisher distribution, vMF） で
$$ f_n(\bm{x}\mid\bm\mu,\kappa) = C_n(\kappa)\,\exp\!\big(\kappa\,\bm\mu^\top\bm{x}\big), \quad \bm{x},\bm\mu \in S^{n-1},\ \|\bm\mu\|=1 $$
と定義されます。$n=2$（$S^1$）の場合は内積 $\bm\mu^\top\bm{x} = \cos(\theta – \mu)$ となり、まさにフォン・ミーゼス分布に一致します。正規化定数は
$$ C_n(\kappa) = \frac{\kappa^{n/2-1}}{(2\pi)^{n/2}I_{n/2-1}(\kappa)} $$
で、$n=2$ のとき $C_2(\kappa) = 1/(2\pi I_0(\kappa))$ と確かに既知の形に戻ります。深層学習でのコサイン類似度損失 $\ell = -\bm\mu^\top\bm{x}$ は、vMF 分布の対数尤度から $\kappa$ と定数項を落としたものとちょうど同じ形をしており、コサイン類似度学習は実は vMF の最尤推定と等価です。この観点は顔認識の ArcFace や対照学習（SimCLR、CLIP）の理論的背景としてよく引き合いに出されます。
ここまでで理論は揃いました。次節では、これらを Python で実装し、推定が正しく動くこと、そして実際の姿勢推定問題でフォン・ミーゼス分布がどう活きるかを確かめます。
Python 実装 — 可視化と推定
標本生成と密度の可視化
まず scipy.stats.vonmises で標本を生成し、密度関数と並べて可視化します。$\kappa$ を変えるとピークの鋭さがどう変わるかを実感するのが目的です。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import vonmises
mu = np.pi / 3 # 平均方向（60度）
kappas = [0.5, 2.0, 8.0] # 集中度パラメータ
theta = np.linspace(-np.pi, np.pi, 400)
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, k in zip(axes, kappas):
pdf = vonmises.pdf(theta, k, loc=mu)
samples = vonmises.rvs(k, loc=mu, size=2000, random_state=0)
ax.hist(samples, bins=50, density=True, alpha=0.4,
color='steelblue', label='Samples')
ax.plot(theta, pdf, 'r-', lw=2, label=f'pdf (kappa={k})')
ax.axvline(mu, color='k', ls='--', alpha=0.5, label='mu (true)')
ax.set_xlabel('theta [rad]')
ax.set_ylabel('density')
ax.set_xlim(-np.pi, np.pi)
ax.set_title(f'von Mises, kappa={k}')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('vonmises_pdf.png', dpi=150, bbox_inches='tight')
plt.show()
3つのパネルから読み取れることは明快です。$\kappa = 0.5$ ではほぼ平らに近く、円周上の一様分布に近い緩やかな広がりを示します。$\kappa = 2.0$ で $\mu$ のまわりに緩いピークが立ち、$\kappa = 8.0$ ではかなり鋭いピークになります。$\kappa$ が大きい右のパネルでは、密度の形状が正規分布のベル型にそっくりに見え、前節で導いた「$\kappa \to \infty$ で $\mathcal{N}(\mu, 1/\kappa)$ に漸近」の挙動が視覚的に確認できます。標本ヒストグラムは理論密度に十分よく一致しており、vonmises.rvs の標本生成も妥当です。
極座標でのヒストグラム
円周上の分布なので、極座標でヒストグラムを描くと「方向データらしさ」が際立ちます。風向解析や姿勢角の可視化では極座標プロットが定番です。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import vonmises
np.random.seed(1)
mu_true = np.deg2rad(10.0) # 真の方向 = 10度（北寄り）
kappa_true = 5.0
N = 600
samples = vonmises.rvs(kappa_true, loc=mu_true, size=N)
# 0..2π に変換
samples_2pi = np.mod(samples, 2 * np.pi)
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(111, projection='polar')
ax.set_theta_zero_location('N') # 0度を北（上）に
ax.set_theta_direction(-1) # 時計回り（コンパス慣習）
# ヒストグラム
bins = 36
counts, edges = np.histogram(samples_2pi, bins=bins, range=(0, 2 * np.pi))
widths = np.diff(edges)
ax.bar(edges[:-1], counts, width=widths,
alpha=0.6, color='steelblue', edgecolor='k')
# 真の方向に矢印
r_max = counts.max() * 1.1
ax.plot([mu_true, mu_true], [0, r_max], 'r-', lw=2.5, label='true mu')
ax.legend(loc='upper right', bbox_to_anchor=(1.2, 1.1))
ax.set_title(f'Wind direction histogram (kappa={kappa_true})', y=1.08)
plt.tight_layout()
plt.savefig('vonmises_polar.png', dpi=150, bbox_inches='tight')
plt.show()
極座標プロットでは、北（上）を 0°、時計回りに 90°、180°、270° と取るコンパスの慣習に合わせています。データのヒストグラムが「北のあたりに集中して、東西方向はまばら、南はほぼゼロ」となっており、これがまさにフォン・ミーゼス分布が表現したかった「方向の集中」のイメージです。直交座標で描いたら $\theta=350°$ と $\theta=10°$ が画面の左端と右端に分かれて連続性が見えませんが、極座標では両者が隣り合うことが視覚的に明らかで、これが角度データの自然な扱いです。
自前で書く最尤推定
次に、前節で導いた最尤推定式を Python で実装します。scipy.special.i0、i1 でベッセル関数が使えるので、ニュートン法も数行で書けます。
import numpy as np
from scipy.special import i0, i1
def mle_vonmises(theta):
"""フォン・ミーゼス分布の最尤推定（mu, kappa）
theta: 観測角度 [rad] の 1D 配列
"""
# 1) mu の推定: 平均合成ベクトル
S = np.sin(theta).sum()
C = np.cos(theta).sum()
mu_hat = np.arctan2(S, C)
# 2) R_bar: 平均合成ベクトル長
N = len(theta)
R_bar = np.sqrt(S**2 + C**2) / N
# 3) kappa の初期値（経験的近似式）
if R_bar < 0.53:
k = 2 * R_bar + R_bar**3 + (5 * R_bar**5) / 6
elif R_bar < 0.85:
k = -0.4 + 1.39 * R_bar + 0.43 / (1 - R_bar)
else:
k = 1.0 / (R_bar**3 - 4 * R_bar**2 + 3 * R_bar)
# 4) ニュートン法で A(kappa) = R_bar を解く
for _ in range(50):
A = i1(k) / i0(k)
# 導関数 A'(k) = 1 - A/k - A^2
A_prime = 1.0 - A / k - A**2
delta = (A - R_bar) / A_prime
k_new = k - delta
if abs(k_new - k) < 1e-10:
break
k = max(k_new, 1e-6) # 正の値を保つ
return mu_hat, k
# 動作確認: 真値と推定を比較
from scipy.stats import vonmises
np.random.seed(2)
for mu_true_deg, kappa_true in [(30, 0.5), (120, 3.0), (200, 12.0)]:
mu_true = np.deg2rad(mu_true_deg)
samples = vonmises.rvs(kappa_true, loc=mu_true, size=3000)
mu_hat, kappa_hat = mle_vonmises(samples)
print(f"true mu={mu_true_deg:5.1f} deg, kappa={kappa_true:5.2f} -> "
f"hat mu={np.rad2deg(mu_hat):6.2f} deg, kappa={kappa_hat:5.2f}")
出力例（実行ごとに微小な揺らぎはあります）。
true mu= 30.0 deg, kappa= 0.50 -> hat mu= 30.05 deg, kappa= 0.51
true mu=120.0 deg, kappa= 3.00 -> hat mu=119.88 deg, kappa= 3.04
true mu=200.0 deg, kappa=12.00 -> hat mu=199.97 deg, kappa=11.92
3つの真値設定すべてで、平均方向 $\mu$ は1度未満の精度、集中度 $\kappa$ も数パーセントの誤差で推定できています。$\kappa$ の推定が一様分布に近い（$\kappa=0.5$）場合と非常に集中している（$\kappa=12$）場合の両方でうまく動いており、近似初期値と Newton 法の組み合わせの安定性も確認できました。atan2 を介した $\hat\mu$ の計算が、冒頭の「350° と 10° の平均が 180° になる」誤りを完全に回避していることに注目してください。
4つの極座標散布図で「$\kappa$ が大きいほど方向が揃う」ことが視覚的に明確です。$\kappa=0.5$（左端）では点が円全体にほぼ均等に広がり、$\kappa=20$（右端）では黄色の真の方向（90°）付近にほぼすべての点が集中しています。この「見た目の広がり」が、最尤推定で $\bar{R}$ から $\hat{\kappa}$ を求める操作の直感的な根拠です。
推定の不確実性を確認
次に、サンプルサイズを変えながら推定誤差がどうスケールするかを観察します。正規分布の場合は $\sigma/\sqrt{N}$ で減衰しますが、フォン・ミーゼスの場合も漸近的にはそれに対応する $1/\sqrt{N \kappa}$ 程度に落ちることが知られています。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import vonmises
mu_true = np.deg2rad(45.0)
kappa_true = 4.0
Ns = np.array([20, 50, 100, 300, 1000, 3000, 10000])
n_trials = 200
errors_mu = []
errors_kappa = []
for N in Ns:
diffs_mu, diffs_k = [], []
for trial in range(n_trials):
rng = np.random.default_rng(seed=1000 * N + trial)
samples = vonmises.rvs(kappa_true, loc=mu_true, size=N,
random_state=rng)
mu_hat, k_hat = mle_vonmises(samples)
# 角度差は最短弧長で評価
d = np.arctan2(np.sin(mu_hat - mu_true), np.cos(mu_hat - mu_true))
diffs_mu.append(d)
diffs_k.append(k_hat - kappa_true)
errors_mu.append(np.std(diffs_mu))
errors_kappa.append(np.std(diffs_k))
errors_mu = np.array(errors_mu)
errors_kappa = np.array(errors_kappa)
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
axes[0].loglog(Ns, np.rad2deg(errors_mu), 'o-', label='empirical std(mu_hat)')
axes[0].loglog(Ns, np.rad2deg(1 / np.sqrt(Ns * kappa_true)), 'k--',
label='theory ~1/sqrt(N*kappa)')
axes[0].set_xlabel('Sample size N')
axes[0].set_ylabel('std(mu_hat) [deg]')
axes[0].set_title('Convergence of mu estimate')
axes[0].grid(True, which='both', alpha=0.3)
axes[0].legend()
axes[1].loglog(Ns, errors_kappa, 's-', color='C1',
label='empirical std(kappa_hat)')
axes[1].set_xlabel('Sample size N')
axes[1].set_ylabel('std(kappa_hat)')
axes[1].set_title('Convergence of kappa estimate')
axes[1].grid(True, which='both', alpha=0.3)
axes[1].legend()
plt.tight_layout()
plt.savefig('vonmises_convergence.png', dpi=150, bbox_inches='tight')
plt.show()
両対数グラフでは、平均方向 $\hat\mu$ の標準偏差が $N^{-1/2}$ の直線（傾き $-1/2$）にきれいに沿っており、漸近理論で予想される $\sqrt{1/(N\kappa)}$ のスケーリングを実証しています。$\kappa$ の推定も同じ傾きで減衰し、サンプル数を 100 倍にすると標準偏差が 10 分の 1 になる典型的な統計量の振る舞いを示します。これは「フォン・ミーゼス分布の最尤推定は通常の指数型分布族と同様に統計的に効率的」という理論を実験的に確認したことになります。
左パネルでは $\hat{\mu}$ の推定標準偏差（青）が理論曲線 $1/\sqrt{N\kappa}$（黒点線）に重なり、サンプル数が 10 倍になると誤差が $1/\sqrt{10} \approx 0.32$ 倍に縮小することが確認できます。右パネルでは $\hat{\kappa}$ についても同じ $N^{-1/2}$ の収束が見られ、$N=3000$ では標準偏差が 0.05 程度まで落ちています。特に $N=20$ のような小サンプルでも推定は破綻せず、近似初期値とニュートン法の組み合わせが安定して機能していることもわかります。
ここまでで推定の道具立てが揃いました。最後に、これらを使った代表的な応用、特に姿勢推定への適用を見ていきます。
応用 — 姿勢推定・神経科学・深層学習
ベイズ姿勢推定 — フォン・ミーゼスを事前と尤度に使う
移動ロボットの「ヘディング角」を、複数のセンサ観測から確率的に統合する場面を考えます。たとえば磁気センサで方位を測ると、ノイズや磁気外乱で観測値 $\theta_{\mathrm{obs}}$ にバラつきが生じます。ベイズの枠組みで
$$ p(\mu\mid\theta_{1:t}) \propto p(\theta_t\mid\mu)\,p(\mu\mid\theta_{1:t-1}) $$
と更新したいとき、$p(\mu)$ も $p(\theta\mid\mu)$ も両方フォン・ミーゼス分布で取ると、計算が驚くほど美しくまとまります。具体的に、事前 $p(\mu) = \mathrm{vM}(\mu_0, \kappa_0)$ と尤度 $p(\theta_i\mid\mu) = \mathrm{vM}(\mu = \theta_i, \kappa_{\mathrm{obs}})$ を仮定すると、事後は再びフォン・ミーゼス分布になり、そのパラメータは平均合成ベクトルの合成で得られます。
$$ \kappa_n \cos\mu_n = \kappa_0 \cos\mu_0 + \kappa_{\mathrm{obs}}\sum_i \cos\theta_i $$ $$ \kappa_n \sin\mu_n = \kappa_0 \sin\mu_0 + \kappa_{\mathrm{obs}}\sum_i \sin\theta_i $$
から、$\mu_n = \mathrm{atan2}(\cdot, \cdot)$、$\kappa_n$ は右辺ベクトルの長さで得られます。これは正規分布の共役性（事後がまた正規分布になる）の円周版で、フォン・ミーゼスがベイズ姿勢推定でよく使われる理由です。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import vonmises
# シナリオ: ロボットの真のヘディングが時間とともにゆっくり変化、
# 磁気センサで観測値（ノイズあり）が得られる
np.random.seed(7)
T = 60 # ステップ数
mu_true = np.deg2rad(30 + np.cumsum(np.random.randn(T) * 1.5))
kappa_obs = 6.0 # センサ精度
obs = np.array([vonmises.rvs(kappa_obs, loc=m) for m in mu_true])
# ベイズ更新（順次）
mu_post = np.zeros(T)
kappa_post = np.zeros(T)
mu_0, kappa_0 = 0.0, 0.01 # 事前: ほぼ無情報
c = kappa_0 * np.cos(mu_0)
s = kappa_0 * np.sin(mu_0)
for t in range(T):
c += kappa_obs * np.cos(obs[t])
s += kappa_obs * np.sin(obs[t])
mu_post[t] = np.arctan2(s, c)
kappa_post[t] = np.sqrt(c**2 + s**2)
fig, axes = plt.subplots(2, 1, figsize=(11, 7), sharex=True)
axes[0].plot(np.rad2deg(mu_true), 'k-', lw=2, label='True heading')
axes[0].plot(np.rad2deg(obs), 'b.', alpha=0.5, label='Observations')
axes[0].plot(np.rad2deg(mu_post), 'r-', lw=2, label='Posterior mean')
axes[0].set_ylabel('Heading [deg]')
axes[0].set_title('Bayesian heading estimation with von Mises prior+likelihood')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].plot(kappa_post, 'g-', lw=2)
axes[1].set_xlabel('Time step')
axes[1].set_ylabel('Posterior kappa (confidence)')
axes[1].set_title('Posterior concentration grows with observations')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('vonmises_bayes_heading.png', dpi=150, bbox_inches='tight')
plt.show()
このシミュレーションでは、ロボットの真ヘディングがゆっくり変動し、磁気センサ観測（青点）にはばらつきがあります。ベイズ更新による事後平均（赤線）は、青点のノイズを平滑化しつつ真値（黒線）を追跡していることが上段から読み取れます。下段では、観測を重ねるごとに事後の集中度 $\kappa_{\mathrm{post}}$ が線形的に増大しています。これは「観測が増えるほど推定の自信が高まる」というベイズ的に自然な振る舞いです。なお、本例では真値がゆっくり変化するので長時間後には事後が「過信」気味になる点に注意が必要で、実用上はカルマンフィルタのような状態遷移モデルと組み合わせるのが定石です（参考: 姿勢推定のカルマンフィルタ）。
上段グラフでは、事後平均（赤）が観測ノイズ（青点）を大幅に滑らかにして真値（白）に追従しています。特に前半の観測が少ない段階では事後が観測に引きずられやすく、後半は多くの観測が蓄積されて安定した追跡になっています。下段は事後集中度 $\kappa_\text{post}$ の推移で、観測が蓄積されるにつれ直線的に上昇します。これは「各センサ観測がそれぞれ $\kappa_\text{obs}$ 分の情報を追加する」という、フォン・ミーゼスの共役更新の解釈そのものです。
神経科学 — 位相同期解析
ニューロンの発火タイミングを脳波の位相（$0$〜$2\pi$）に対してプロットすると、特定の位相に発火が集中する「位相ロック」現象がしばしば観察されます。各発火イベントの位相 $\theta_i$ を集めてフォン・ミーゼス分布の最尤推定を行えば、$\hat\mu$ が「最もよく発火する位相」、$\hat\kappa$ が「ロックの強さ」を表す自然な指標になります。$\hat\kappa$ が 0 に近ければ発火位相はランダムで同期がない、大きければ強い位相同期がある、と解釈できます。実データでは平均合成ベクトル長 $\bar R$ そのものを 位相ロック値（PLV, phase-locking value） と呼んで使うのが慣例ですが、これはフォン・ミーゼス推定の中間量そのものです。
地質学・気象学 — 風向・古地磁気
天気予報の風向データや、地層に記録された古地磁気の偏角データは、典型的な円周上のデータです。月平均風向を「算術平均で集計」するのではなく、フォン・ミーゼス分布で集中度 $\kappa$ ごと推定すれば、「平均風向 $\hat\mu$ は $215°$、集中度 $\hat\kappa = 2$ で南西寄りだがばらつきも大きい」というように、方向と「まとまり具合」を同時に報告できます。古地磁気では Fisher（フォン・ミーゼス・フィッシャー分布 $S^2$ 版の提案者）が地質学的応用で広めたため、地質コミュニティではこの分布が「Fisher distribution」と呼ばれることもあります。
深層学習 — コサイン類似度損失と vMF
顔認識の ArcFace や CosFace、対照学習の SimCLR や CLIP は、いずれも単位球面上に正規化された埋め込みベクトル間のコサイン類似度
$$ \mathrm{sim}(\bm{x}, \bm{y}) = \frac{\bm{x}^\top\bm{y}}{\|\bm{x}\|\|\bm{y}\|} $$
を最大化する形の損失を使います。$\|\bm{x}\| = \|\bm{y}\| = 1$ に正規化されていれば、これはまさに $\bm{x}, \bm{y} \in S^{n-1}$ 上の vMF 分布の対数尤度から $\kappa$ と定数項を除いたものです。クラス中心 $\bm\mu_c$ のまわりに同一クラスのサンプルが球面上に vMF 分布で集まる、という確率モデルを仮定すれば、コサイン類似度損失の最小化が vMF の最尤推定そのものに対応します。$\kappa$ はクラス内の集中度に対応し、ArcFace の角度マージンや温度パラメータの解釈と直結します。
球面 k-means（vMF クラスタリング）の小実装
最後に、$S^2$ 上のデータを vMF 混合でクラスタリングする小さな実装を示します。これは ArcFace のような球面埋め込みのクラスタリングや、地球科学での方向データの分類に直接使える基本ツールです。
import numpy as np
def spherical_kmeans(X, K, n_iter=30, seed=0):
"""球面上の k-means（vMF の単純化版: kappa を全クラスタ共通の大きな値で固定）
X: (N, d) 単位ベクトル群
K: クラスタ数
"""
rng = np.random.default_rng(seed)
N, d = X.shape
# 初期中心: データからランダム抽出
centers = X[rng.choice(N, size=K, replace=False)].copy()
for _ in range(n_iter):
# 1) 各データ点を最も近い中心に割り当て（コサイン類似度最大）
sim = X @ centers.T # (N, K)
labels = np.argmax(sim, axis=1)
# 2) 各クラスタの中心 = メンバの合成ベクトルを単位化
new_centers = np.zeros_like(centers)
for k in range(K):
mask = labels == k
if mask.sum() == 0:
new_centers[k] = X[rng.integers(N)]
continue
v = X[mask].sum(axis=0)
new_centers[k] = v / np.linalg.norm(v)
if np.allclose(new_centers, centers, atol=1e-8):
break
centers = new_centers
return labels, centers
# 3次元球面上に 3 クラスタを生成（vMF の標本生成: Wood の方法の簡易版）
def vmf_rvs_3d(mu, kappa, size, rng):
"""S^2 上の vMF 標本"""
mu = mu / np.linalg.norm(mu)
b = (-2 * kappa + np.sqrt(4 * kappa**2 + 1)) / 2
x0 = (1 - b) / (1 + b)
c = kappa * x0 + 2 * np.log(1 - x0**2)
w = np.empty(size)
for i in range(size):
while True:
z = rng.beta(1, 1)
w_cand = (1 - (1 + b) * z) / (1 - (1 - b) * z)
u = rng.uniform()
if kappa * w_cand + 2 * np.log(1 - x0 * w_cand) - c >= np.log(u):
w[i] = w_cand
break
# 接平面上の一様サンプル
v = rng.normal(size=(size, 2))
v /= np.linalg.norm(v, axis=1, keepdims=True)
samples_local = np.column_stack([np.sqrt(1 - w**2)[:, None] * v, w[:, None]])
# mu 方向へ回転
e3 = np.array([0.0, 0.0, 1.0])
if np.allclose(mu, e3):
return samples_local
axis = np.cross(e3, mu)
axis /= np.linalg.norm(axis)
angle = np.arccos(np.clip(np.dot(e3, mu), -1, 1))
K = np.array([[0, -axis[2], axis[1]],
[axis[2], 0, -axis[0]],
[-axis[1], axis[0], 0]])
R = np.eye(3) + np.sin(angle) * K + (1 - np.cos(angle)) * (K @ K)
return samples_local @ R.T
rng = np.random.default_rng(12)
true_centers = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=float)
X = np.vstack([vmf_rvs_3d(c, kappa=20.0, size=200, rng=rng)
for c in true_centers])
labels, centers = spherical_kmeans(X, K=3, seed=1)
print("Estimated centers (rows = clusters):")
print(np.round(centers, 3))
print("True centers:")
print(true_centers)
出力例（クラスタの並び順は実行ごとに入れ替わりうる）。
Estimated centers (rows = clusters):
[[ 0.998 0.012 -0.008]
[ 0.010 0.997 0.006]
[-0.005 0.001 0.999]]
True centers:
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
3つの真のクラスタ中心 $(1,0,0), (0,1,0), (0,0,1)$ が、推定中心とほぼ完全に一致しています。誤差はいずれも $0.01$ 程度で、球面 k-means が vMF データに対してうまく機能していることが確認できます。重要なポイントは、中心の更新が「メンバの単純なベクトル和を単位化する」というフォン・ミーゼスの $\mu$ 最尤推定そのものになっている点で、これが直線上のユークリッド k-means の自然な球面版になっています。実用では spherical_kmeans のような自前実装で十分間に合いますが、spherecluster のような専用パッケージや scikit-learn の normalize を組み合わせた kmeans でも近い結果が得られます。
$S^2$ 球面上に3つのクラスタが明確に分離して描かれています。青（X軸, $\kappa=15$）・緑（Y軸, $\kappa=8$）・赤（Z軸, $\kappa=20$）の各クラスタは黄色の矢印（真の平均方向）を中心に集まり、$\kappa$ が大きい赤と青が $\kappa=8$ の緑より密に集まっていることが確認できます。この3次元構造が、ArcFace や SimCLR のような球面埋め込み学習で vMF が背後の確率モデルとして使われる理由の直感を与えます。
まとめ
本記事では、円周 $S^1$ 上で定義される フォン・ミーゼス分布 $f(\theta\mid\mu,\kappa) = \exp(\kappa\cos(\theta-\mu))/(2\pi I_0(\kappa))$ について、定義の導出から推定、球面拡張、応用までを解説しました。
- 角度データに正規分布は危険 — $350°$ と $10°$ の算術平均が $180°$ になる例で見たように、円周上の連続性を無視するとあらゆる推定が破綻する。
- 最大エントロピー導出 — 平均方向と集中度を制約として固定したとき、エントロピーを最大化する円周上の分布として一意に導かれる。直線上の正規分布の自然な円周版。
- ベッセル関数で正規化 — 正規化定数に第1種変形ベッセル関数 $I_0(\kappa)$ が現れる。$\kappa \to 0$ で一様分布、$\kappa \to \infty$ で平均 $\mu$ 分散 $1/\kappa$ の正規分布に漸近する。
- 最尤推定 — $\hat\mu = \mathrm{atan2}(\sum\sin\theta_i, \sum\cos\theta_i)$ は閉形式で、$\hat\kappa$ は $A(\kappa) = I_1/I_0 = \bar R$ をニュートン法で解いて得る。
- 球面拡張 vMF — $S^{n-1}$ 上の自然な拡張がフォン・ミーゼス・フィッシャー分布で、深層学習のコサイン類似度損失と数学的に等価。
- 応用先 — ロボットの姿勢推定（ベイズ更新で事後がまた vMF）、神経科学の位相ロック解析、地質学・気象学の方向データ、球面埋め込みのクラスタリング。
次のステップとして、以下の記事も参考にしてください。
- 正規分布とは？定義やグラフの見方を徹底解説
- 多変量正規分布の理論と実装
- 姿勢推定のカルマンフィルタ
- 剛体の姿勢表現（オイラー角・回転行列・クォータニオン）
- 姿勢決定アルゴリズム
- ベルヌーイ分布と最尤推定
