-
Notifications
You must be signed in to change notification settings - Fork 1
/
circ_vmum_est_mm.m
80 lines (60 loc) · 2.82 KB
/
circ_vmum_est_mm.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
function [mu_hat, k_hat, p1_hat, p2_hat, p3_hat] = circ_vmum_est_mm(data, mu)
%CIRC_VMUM_EST_MM calculates MM estimate of vMUM model parameters
%
% Audio Circular Statistics (ACS) library
% Copyright 2016 Enzo De Sena
%% Asserts
assert(iscolumn(data));
assert(length(data) >= 1);
phi = mod(data*2, 2*pi);
%% Estimate mu
if nargin == 2 % if mu is known, then mu_hat = mu
assert(isscalar(mu));
mu_hat = mu;
else
mu_hat = circ_mean(phi)/2;
end
%% Estimate k
ac2w = mean(cos(2*(phi-2*mu_hat)));
ac1w = mean(cos(1*(phi-2*mu_hat)));
%options = optimoptions('fsolve','TolFun',1E-9);
options = optimset('Display', 'off');
k_hat = fsolve(@(k)besseli(4,k)./besseli(2,k)*ac1w-ac2w, 1, options);
%% Assumes prior knowledge that k<100
%k_hat = fsolve(@(k)besseli(4,k)./besseli(2,k)-min(abs(ac2w/ac1w), 0.94), 1, options);
%% Works worse because rho2w>ac2w
% bc2w = mean(sin(2*(phi-2*mu_hat))); rho2w = sqrt(ac2w.^2+bc2w.^2);
% k_hat = fsolve(@(k)ac1w.*besseli(4,k)./besseli(2,k)-rho2w, 1, options);
%k_hat = fsolve(@(k)ac1w*besseli(4,k)./besseli(2,k)-ac2w, 1, options);
%% Version that aims to remove bias
%k_hat = fsolve(@(k)ac1w*besseli(4,k)./besseli(2,k)-ac2w, 1)-4.54E-5*exp(15.06*abs(ac2w/ac1w));
%% Doesn't work because besseli(4,k) and besseli(2,k) are very large numbers
%k_hat = fsolve(@(k)abs(ac1w)*besseli(4,k)-besseli(2,k)*abs(ac2w), 1);
p_ol = min(abs(ac1w)*besseli(0,k_hat)/besseli(2,k_hat), 1);
%% Sequential correction
k_hat = fsolve(@(k)besseli(2,k)./besseli(0,k)*p_ol-(abs(ac1w)), 1, options);
p_ol = min(abs(ac1w)*besseli(0,k_hat)/besseli(2,k_hat), 1);
k_hat = fsolve(@(k)besseli(2,k)./besseli(0,k)*p_ol-(abs(ac1w)), 1, options);
p_ol = min(abs(ac1w)*besseli(0,k_hat)/besseli(2,k_hat), 1);
k_hat = fsolve(@(k)besseli(2,k)./besseli(0,k)*p_ol-(abs(ac1w)), 1, options);
%% Works worse than other methods
% [mean(cos(1*(phi-2*mu_hat))),mean(cos(2*(phi-2*mu_hat))),mean(cos(3*(phi-2*mu_hat))),mean(cos(4*(phi-2*mu_hat)))]
% [mean(sin(1*(phi-2*mu_hat))),mean(sin(2*(phi-2*mu_hat))),mean(sin(3*(phi-2*mu_hat))),mean(sin(4*(phi-2*mu_hat)))]
% params = fmincon(@(params)(...
% (besseli(2,params(1))./besseli(0,params(1))*params(2)-mean(cos(1*(phi-2*mu_hat)))).^2+...
% (besseli(4,params(1))./besseli(0,params(1))*params(2)-mean(cos(2*(phi-2*mu_hat)))).^2+...
% (besseli(6,params(1))./besseli(0,params(1))*params(2)-mean(cos(3*(phi-2*mu_hat)))).^2+...
% (besseli(8,params(1))./besseli(0,params(1))*params(2)-mean(cos(4*(phi-2*mu_hat)))).^2), [1, 0.5],...
% [],[],[],[],[0, 0],[inf,1]);
%
% k_hat = params(1);
% p_ol = params(2);
%% Estimate p1 and p2
a1 = mean(cos(data-mu_hat));
p1_hat = min(abs(p_ol + a1*besseli(0,k_hat)/besseli(1,k_hat))/2, p_ol);
p2_hat = p_ol - p1_hat;
p3_hat = 1-p1_hat-p2_hat;
%% Asserts
assert(p1_hat>=0 && p1_hat<=1);
assert(p2_hat>=0 && p2_hat<=1);
assert((1-p1_hat-p2_hat)>=0 && (1-p1_hat-p2_hat)<=1);