-
Notifications
You must be signed in to change notification settings - Fork 2
/
EstimateParametersHMRFMCEM.m
46 lines (43 loc) · 1.5 KB
/
EstimateParametersHMRFMCEM.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
%% Îöåíèâàåò ïàðàìåòðû beta, mu, kappa ïî âûáîðêå èç ñõåìû Ãèááñà
%---input---------------------------------------------------------
% X: 4-õ ìåðíàÿ ìàòðèöà èñõîäíûõ äàííûõ
% Y: ìàòðèöà ðàçìåðíîñòè MxN
% k: êîëè÷åñòâî êëàññîâ ñåãìåíòàöèè
% p: ðàçìåðíîñòü äàííûõ
% beta: òåêóùåå çíà÷åíèå ïàðàìåòðà ìîäåëè Ïîòòñà
% mu: òåêóùàÿ ìàòðèöà ïàðàìåòðîâ äëÿ vMF
% kappa: òåêóùèé âåêòîð ïàðàìåòðîâ äëÿ vMF
%---output--------------------------------------------------------
% beta: íîâîå çíà÷åíèå ïàðàìåòðà, ñêàëÿð
% mu: íîâûå çíà÷åíèÿ ïàðàìåòðîâ, ìàòðèöà LxP
% kappa: íîâûå çíà÷åíèÿ ïàðàìåòðîâ, âåêòîð 1xL
function [beta, mu, kappa] = EstimateParametersHMRFMCEM(X, Y, k, p, beta, mu, kappa)
M = size(Y, 1);
for l=1:k
R = zeros(1, p);
N = 0;
for i=1:M
R = R + sum(X(Y(i, :)==l, :));
N = N + sum(Y(i, :)==l);
end
R_mod = sqrt(sum(R .^ 2));
if R_mod ~= 0
mu(l, :) = R / R_mod;
end
eqn = @(x) besseli(p/2, x, 1)/besseli(p/2-1, x, 1) * N / M - R_mod / M;
initval = kappa(l);
if isnan(eqn(initval))
fprintf('WARNING: function is NaN in initial point\n')
initval = 10;
end
opts1 = optimset('display','off');
kappa(l) = lsqnonlin(eqn, initval, 0, Inf, opts1);
if kappa(l) < 0
kappa(l) = abs(kappa(l));
fprintf('WARNING: Lower than zero solution was found\n')
end
if isnan(kappa(l))
kappa(l) = initval;
fprintf('WARNING: NaN solution was found\n')
end
end