von Mises distribution(馮·米賽斯分布)的隨機模擬
與參數估計的筆記（一）
1.von Mises distribution概述
在概率論和方向統計中，von-Mises 分布（也稱為圓形正態分布或 Tikhonov 分布）是圓上的連續概率分布。 它是環繞正態分布的近似值，它是正態分布的循環模擬。 圓上的自由擴散角 θ 是一個包裹的正態分布隨機變量，其展開方差隨時間線性增長。 另一方面，von-Mises 分布是諧波勢中圓上漂移和擴散過程的平穩分布，即具有優選方向。 von-Mises 分布是指定第一個圓形矩的實部和虛部時圓形數據的最大熵分布。 von-Mises 分布是 von Mises-Fisher 分布在 N 維球面上的一個特例（如圖1所示）。
圖1 von Mises Distribution 與 Fisher Distribution 的關系圖
1.1 von Mises分布的定義
對於角度x的von Mises分布的概率密度PDF由下式給出：
\[f\left( {x\left| {\mu ,\kappa } \right.} \right) = \frac{{{e^{\kappa \cos \left( {x - \mu } \right)}}}}{{2\pi {I_0}\left( \kappa \right)}} \]
\({{I_0}\left( \kappa \right)}\) 為零階修正貝塞爾公式。其中參數 $\mu $ 和\(\kappa\)可以近似理解為正態分布的\({{N}\left( \mu ,\delta \right)}\)中的均值\(\mu\)和標准差\(\delta\)
-
\(\mu\) 其中是衡量方位角\(x\)的平均值，整個分布的平均體現；
-
\(\kappa\) 其中是衡量方位角\({x}\)的集中程度\(\left(\kappa>0 \right)\)
當 \(\kappa\)趨近於零，分布足夠分散，不足以滿足角度\(x\)均值的存在性，方位角方向滿足，此時的PDF滿足均勻分布，及
\[f\left( x \right) = \frac{1}{{2\pi }} \]
當\(\kappa\)趨近於無窮大時，方位角的集中程度很高，其近似於均值為\(\mu\), 方差為\(\frac{1}{\kappa}\)的高斯分布：
\[f\left( x \right) \approx \frac{{\sqrt \kappa }}{{\sqrt {2\pi } }}\exp \left( { - \frac{{\kappa {{\left( {x - \mu } \right)}^2}}}{2}} \right) \]
其von-Mises分布的PDF可以寫成如下貝塞爾函數的級數序列：
\[f(x\mid \mu ,\kappa )={\frac {1}{2\pi }}\left(1+{\frac {2}{I_{0}(\kappa )}}\sum _{j=1}^{\infty }I_{j}(\kappa )\cos[j(x-\mu )]\right) \]
其中 \(I_{j}(\kappa )\)是\(j\)階改進的貝塞爾函數。
圖 2 von mises分布的PDF圖
其中CDF的公式我們無法獲取解析解，但可以對其Bessel函數展開形式進行積分得到其數值解得到如下式表達：
\[\Phi (x\mid \mu ,\kappa )=\int f(t\mid \mu ,\kappa )\,dt={\frac {1}{2\pi }}\left(x+{\frac {2}{I_{0}(\kappa )}}\sum _{j=1}^{\infty }I_{j}(\kappa ){\frac {\sin[j(x-\mu )]}{j}}\right). \]
其CDF分布如下圖表示：
圖 3 von mises分布的CDF圖
其中的von Mises分布的統計特征參數如下表所示：
表 一 von Mises分布統計參數表
1.2 von Mises 分布的隨機模擬代碼（matlab）
其中模擬與估計代碼來源於Philipp Berens 和 Marc J. Velasco的Circular Statistics Toolbox
function alpha = circ_vmrnd(theta, kappa, n)
% alpha = circ_vmrnd(theta, kappa, n)
% Simulates n random angles from a von Mises distribution, with preferred
% direction thetahat and concentration parameter kappa.
%
% Input:
% [theta preferred direction, default is 0]
% [kappa width, default is 1]
% [n number of samples, default is 10]
%
% If n is a vector with two entries (e.g. [2 10]), the function creates
% a matrix output with the respective dimensionality.
%
% Output:
% alpha samples from von Mises distribution
%
%
% References:
% Statistical analysis of circular data, Fisher, sec. 3.3.6, p. 49
%
% Circular Statistics Toolbox for Matlab
% By Philipp Berens and Marc J. Velasco, 2009
% velasco@ccs.fau.edu
% default parameter
if nargin < 3
n = 10;
end
if nargin < 2
kappa = 1;
end
if nargin < 1
theta = 0;
end
if numel(n) > 2
error('n must be a scalar or two-entry vector!')
elseif numel(n) == 2
m = n;
n = n(1) * n(2);
end
% if kappa is small, treat as uniform distribution
if kappa < 1e-6
alpha = 2*pi*rand(n,1);
return
end
% other cases
a = 1 + sqrt((1+4*kappa.^2));
b = (a - sqrt(2*a))/(2*kappa);
r = (1 + b^2)/(2*b);
alpha = zeros(n,1);
for j = 1:n
while true
u = rand(3,1);
z = cos(pi*u(1));
f = (1+r*z)/(r+z);
c = kappa*(r-f);
if u(2) < c * (2-c) || ~(log(c)-log(u(2)) + 1 -c < 0)
break
end
end
alpha(j) = theta + sign(u(3) - 0.5) * acos(f);
alpha(j) = angle(exp(1i*alpha(j)));
end
if exist('m','var')
alpha = reshape(alpha,m(1),m(2));
end
1.3 von Mises 分布的參數估計代碼（matlab）
function [thetahat kappa] = circ_vmpar(alpha,w,d)
% r = circ_vmpar(alpha, w, d)
% Estimate the parameters of a von Mises distribution.
%
% Input:
% alpha sample of angles in radians
% [w number of incidences in case of binned angle data]
% [d spacing of bin centers for binned data, if supplied
% correction factor is used to correct for bias in
% estimation of r, in radians (!)]
%
% Output:
% thetahat preferred direction
% kappa concentration parameter
%
% PHB 3/23/2009
%
% References:
% Statistical analysis of circular data, N.I. Fisher
%
% Circular Statistics Toolbox for Matlab
% By Philipp Berens, 2009
% berens@tuebingen.mpg.de
alpha = alpha(:);
if nargin < 2
w = ones(size(alpha));
end
if nargin < 3
d = 0;
end
r = circ_r(alpha,w,d);
kappa = circ_kappa(r);
thetahat = circ_mean(alpha,w);
其中circ_r,circ_kappa和circ_mean的子代碼如下：
function r = circ_r(alpha, w, d, dim)
% r = circ_r(alpha, w, d)
% Computes mean resultant vector length for circular data.
%
% Input:
% alpha sample of angles in radians
% [w number of incidences in case of binned angle data]
% [d spacing of bin centers for binned data, if supplied
% correction factor is used to correct for bias in
% estimation of r, in radians (!)]
% [dim compute along this dimension, default is 1]
%
% If dim argument is specified, all other optional arguments can be
% left empty: circ_r(alpha, [], [], dim)
%
% Output:
% r mean resultant length
%
% PHB 7/6/2008
%
% References:
% Statistical analysis of circular data, N.I. Fisher
% Topics in circular statistics, S.R. Jammalamadaka et al.
% Biostatistical Analysis, J. H. Zar
%
% Circular Statistics Toolbox for Matlab
% By Philipp Berens, 2009
% berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
if nargin < 4
dim = 1;
end
if nargin < 2 || isempty(w)
% if no specific weighting has been specified
% assume no binning has taken place
w = ones(size(alpha));
else
if size(w,2) ~= size(alpha,2) || size(w,1) ~= size(alpha,1)
error('Input dimensions do not match');
end
end
if nargin < 3 || isempty(d)
% per default do not apply correct for binned data
d = 0;
end
% compute weighted sum of cos and sin of angles
r = sum(w.*exp(1i*alpha),dim);
% obtain length
r = abs(r)./sum(w,dim);
% for data with known spacing, apply correction factor to correct for bias
% in the estimation of r (see Zar, p. 601, equ. 26.16)
if d ~= 0
c = d/2/sin(d/2);
r = c*r;
end
function kappa = circ_kappa(alpha,w)
%
% kappa = circ_kappa(alpha,[w])
% Computes an approximation to the ML estimate of the concentration
% parameter kappa of the von Mises distribution.
%
% Input:
% alpha angles in radians OR alpha is length resultant
% [w number of incidences in case of binned angle data]
%
% Output:
% kappa estimated value of kappa
%
% References:
% Statistical analysis of circular data, Fisher, equation p. 88
%
% Circular Statistics Toolbox for Matlab
% By Philipp Berens, 2009
% berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
alpha = alpha(:);
if nargin<2
% if no specific weighting has been specified
% assume no binning has taken place
w = ones(size(alpha));
else
if size(w,2) > size(w,1)
w = w';
end
end
N = length(alpha);
if N>1
R = circ_r(alpha,w);
else
R = alpha;
end
if R < 0.53
kappa = 2*R + R^3 + 5*R^5/6;
elseif R>=0.53 && R<0.85
kappa = -.4 + 1.39*R + 0.43/(1-R);
else
kappa = 1/(R^3 - 4*R^2 + 3*R);
end
if N<15 && N>1
if kappa < 2
kappa = max(kappa-2*(N*kappa)^-1,0);
else
kappa = (N-1)^3*kappa/(N^3+N);
end
end
function [mu ul ll] = circ_mean(alpha, w, dim)
%
% mu = circ_mean(alpha, w)
% Computes the mean direction for circular data.
%
% Input:
% alpha sample of angles in radians
% [w weightings in case of binned angle data]
% [dim compute along this dimension, default is 1]
%
% If dim argument is specified, all other optional arguments can be
% left empty: circ_mean(alpha, [], dim)
%
% Output:
% mu mean direction
% ul upper 95% confidence limit
% ll lower 95% confidence limit
%
% PHB 7/6/2008
%
% References:
% Statistical analysis of circular data, N. I. Fisher
% Topics in circular statistics, S. R. Jammalamadaka et al.
% Biostatistical Analysis, J. H. Zar
%
% Circular Statistics Toolbox for Matlab
% By Philipp Berens, 2009
% berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
if nargin < 3
dim = 1;
end
if nargin < 2 || isempty(w)
% if no specific weighting has been specified
% assume no binning has taken place
w = ones(size(alpha));
else
if size(w,2) ~= size(alpha,2) || size(w,1) ~= size(alpha,1)
error('Input dimensions do not match');
end
end
% compute weighted sum of cos and sin of angles
r = sum(w.*exp(1i*alpha),dim);
% obtain mean by
mu = angle(r);
% confidence limits if desired
if nargout > 1
t = circ_confmean(alpha,0.05,w,[],dim);
ul = mu + t;
ll = mu - t;
end
function t = circ_confmean(alpha, xi, w, d, dim)
%
% t = circ_mean(alpha, xi, w, d, dim)
% Computes the confidence limits on the mean for circular data.
%
% Input:
% alpha sample of angles in radians
% [xi (1-xi)-confidence limits are computed, default 0.05]
% [w number of incidences in case of binned angle data]
% [d spacing of bin centers for binned data, if supplied
% correction factor is used to correct for bias in
% estimation of r, in radians (!)]
% [dim compute along this dimension, default is 1]
%
% Output:
% t mean +- d yields upper/lower (1-xi)% confidence limit
%
% PHB 7/6/2008
%
% References:
% Statistical analysis of circular data, N. I. Fisher
% Topics in circular statistics, S. R. Jammalamadaka et al.
% Biostatistical Analysis, J. H. Zar
%
% Circular Statistics Toolbox for Matlab
% By Philipp Berens, 2009
% berens@tuebingen.mpg.de - www.kyb.mpg.de/~berens/circStat.html
if nargin < 5
dim = 1;
end
if nargin < 4 || isempty(d)
% per default do not apply correct for binned data
d = 0;
end
if nargin < 3 || isempty(w)
% if no specific weighting has been specified
% assume no binning has taken place
w = ones(size(alpha));
else
if size(w,2) ~= size(alpha,2) || size(w,1) ~= size(alpha,1)
error('Input dimensions do not match');
end
end
% set confidence limit size to default
if nargin < 2 || isempty(xi)
xi = 0.05;
end
% compute ingredients for conf. lim.
r = circ_r(alpha,w,d,dim);
n = sum(w,dim);
R = n.*r;
c2 = chi2inv((1-xi),1);
% check for resultant vector length and select appropriate formula
t = zeros(size(r));
for i = 1:numel(r)
if r(i) < .9 && r(i) > sqrt(c2/2/n(i))
t(i) = sqrt((2*n(i)*(2*R(i)^2-n(i)*c2))/(4*n(i)-c2)); % equ. 26.24
elseif r(i) >= .9
t(i) = sqrt(n(i)^2-(n(i)^2-R(i)^2)*exp(c2/n(i))); % equ. 26.25
else
t(i) = NaN;
warning('Requirements for confidence levels not met.');
end
end
% apply final transform
t = acos(t./R);
1.4 參考引用
（1）https://en.wikipedia.org/wiki/Von_Mises%E2%80%93Fisher_distribution
(2) https://blog.csdn.net/weixin_42576437/article/details/107124117