フォンミーゼス分布サンプリング

確率統計
フォンミーゼス分布サンプリング; von Mises Distribution Sampling

概要

フォンミーゼス分布とは以下の式で表現される確率分布。円筒状の正規分布と表現されることが多い。余談であるが球面上の正規分布(フォンミーゼスフィッシャー分布)もある。

$$f(\theta)=\frac{exp(\kappa cos(\theta-\mu))}{2 \pi I_0(\kappa)} $$

フォンミーゼス分布に従う確率変数の生成にはBarabesi(1995)[1]のアルゴリズムを用いる。単純に正規分布を用いて生成出来ないので注意が必要。

ソースコード


namespace ExRandom.Continuous {
    public class VonMisesRandom : Random{
        readonly MT19937 mt;
        readonly double kappa, mu, s;
        
        public VonMisesRandom(MT19937 mt, double kappa, double mu) {
            if(mt == null) {
                throw new ArgumentNullException();
            }

            this.mt = mt;
            this.kappa = kappa;
            this.mu = RoundMu(mu);
            this.s = (kappa > 1.3) ? (1 / Math.Sqrt(kappa)) : (Math.PI * Math.Exp(-kappa));
        } 

        public override double Next() {
            double t;

            while(true) {
                double r1 = mt.NextDouble_OpenInterval01(), r2 = mt.NextDouble_OpenInterval01();
                t = s * (2 * r2 - 1) / r1;

                if(Math.Abs(t) > Math.PI) {
                    continue;
                }

                if((kappa * t * t < 4 - 4 * r1) || (kappa * Math.Cos(t) >= 2 * Math.Log(r1) + kappa)) {
                    break;
                }
            }

            return RoundTheta(t + mu);
        }

        private static double RoundTheta(double theta) {
            if(theta > Math.PI) {
                theta -= 2 * Math.PI;
            }
            else if(theta < -Math.PI) {
                theta += 2 * Math.PI;
            }

            return theta;
        }

        private static double RoundMu(double mu) {
            mu += Math.PI;

            if(mu >= 0) {
                mu %= 2 * Math.PI;
            }
            else {
                mu = 2 * Math.PI - ((-mu) % (2 * Math.PI));
            }

            mu -= Math.PI;

            return mu;
        }
    }
}

[1] “Generating Von Mises Variates By The Ratio-of-Uniforms Method”, Lucio Barabesi, Statistica Applicata, vol.7, no.4, 1995