von Mises Distribution (馮·米賽斯分佈)的隨機模擬與引數估計的筆記（二）

1.引數估計運算元分析

​ 在上一節中，我們討論了von Mises Distribution的概率分佈函式PDF和累計分佈函式CDF，並給出了von Mises Distribution的隨機模擬和引數估計matlab程式，其中在此我們就引數估計的細節進行補充。其基於最大似然引數估計運算元，如下表：

​ 來源於《Statistical Distributions》

利用如下改進貝塞爾函式的關係求解引數\(\kappa\)

1.1 \(\mu\)

引數估計分析matlab程式碼

function mu=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
 
%
% 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
% [email protected] - 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);

1.2 \(\kappa\)引數估計的matlab程式碼

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
% [email protected] - 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 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
% [email protected] - 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

2 程式碼效果分析

clc
clear all
close all

theta=pi/2;  %設定模擬引數
kappa=50;
n=3000;

alpha=circ_vmrnd(theta,kappa,n);  %生成制定引數的von-Mises分佈的隨機數

[thetahat1 kappa1]=circ_vmpar(alpha); %對其進行分佈引數進行估計分析

 %繪製模擬資料直方圖
figure(1) 
hist(alpha,100);
xlabel('Angle(弧度)');
ylabel('Frequency');

X = categorical({'Really value','Estimate value'});

 %估計引數與模型引數對比
figure(2)   
subplot(1,2,1)
bar(X,[theta,thetahat1]);
ylabel('theta');

subplot(1,2,2)
bar(X,[kappa,kappa1]);
ylabel('kappa');