フォンミーゼスフィッシャー分布サンプリング

確率統計
フォンミーゼスフィッシャー分布サンプリング; von Mises-Fisher Distribution Sampling

概要

フォンミーゼスフィッシャー分布とは以下の式で表現される確率分布。球面上の正規分布ともいえる方位確率分布の一つ。フォンミーゼスフィッシャー分布は3以上の次元を扱えるが、この記事では3次元のときを扱う。
ここで$\boldsymbol{\mu} $は確率密度が最も密となる方位、$ {\kappa} \gt 0 $は密度を示しこの値が大きいときは$\boldsymbol{\mu} $周辺の確率密度が大きくなり、限りなく0に近いときは球面一様分布になる。

$$ f(\boldsymbol{v}) = \frac{\kappa}{4 \pi sinh(\kappa)} exp(\kappa \boldsymbol{\mu}^T \boldsymbol{v} ), \quad |\boldsymbol{v}|=1 $$

フォンミーゼスフィッシャー分布に従う確率変数の生成は以下の手順で行う。
1. 密度$ {\kappa} $で$\boldsymbol{\mu_0} = (1, 0, 0)^T$に対して密になるように確率変数$\boldsymbol{v_0}$を生成
2. $\boldsymbol{\mu_0}$を$\boldsymbol{\mu} $に回転する四元数を確率変数$\boldsymbol{v_0}$に適用し$\boldsymbol{v}$を得る

ソースコード


namespace ExRandom.MultiVariate {
    public class VonMisesFisherRandom : Random<double>{
        readonly MT19937 mt;
        readonly double kappa, inv_kappa, exp_m2kappa;
        readonly double qri, qrj, qij, qxx, qyy, qzz;

        public VonMisesFisherRandom(MT19937 mt, Vector<double> mu, double kappa = 1) {
            if(mt == null) {
                throw new ArgumentNullException();
            }

            if(mu.Dim != 3) { 
                throw new ArgumentException();
            }
            
            if(!(kappa >= 0)) {
                throw new ArgumentException();
            }

            double inv_norm = 1 / Math.Sqrt(mu.X * mu.X + mu.Y * mu.Y + mu.Z * mu.Z);

            if(double.IsNaN(inv_norm) || double.IsInfinity(inv_norm)) {
                throw new ArgumentException();
            }

            this.mt = mt;
            this.kappa = kappa;
            this.inv_kappa = 1 / kappa;
            this.exp_m2kappa = Math.Exp(-2 * kappa);

            double mx, my, mz, qr, qi, qj, angle, c, s, norm;

            mx = mu.X * inv_norm;
            my = mu.Y * inv_norm;
            mz = mu.Z * inv_norm;

            angle = Math.Acos(mz);
            c = Math.Cos(angle * 0.5);
            s = Math.Sin(angle * 0.5);
            norm = Math.Sqrt(mx * mx + my * my);

            if(norm > 0) {
                qr = c;
                qi = s * -my / norm;
                qj = s * mx / norm;
            }
            else {
                qr = 1;
                qi = 0;
                qj = 0;
            }

            this.qri = qr * qi;
            this.qij = qi * qj; 
            this.qrj = qr * qj;
            this.qxx = qr * qr + qi * qi - qj * qj;
            this.qyy = qr * qr - qi * qi + qj * qj;
            this.qzz = qr * qr - qi * qi - qj * qj;
        }

        public override Vector<double> Next() {
            double r, w, s, theta, x, y, z;

            r = mt.NextDouble();
            w = (kappa > 0) ? (1 + Math.Log(r + (1 - r) * exp_m2kappa) * inv_kappa) : (2 * r - 1);
            s = Math.Sqrt(1 - w * w);
            theta = 2 * Math.PI * mt.NextDouble_OpenInterval1();

            x = s * Math.Cos(theta);
            y = s * Math.Sin(theta);
            z = w;

            return new Vector<double>(x * qxx + 2 * (y * qij + z * qrj), y * qyy - 2 * (z * qri - x * qij), z * qzz - 2 * (x * qrj - y * qri));
        }
    }
}

関連項目

メルセンヌ・ツイスタ
各種確率分布サンプリング基本クラス