-
Notifications
You must be signed in to change notification settings - Fork 2
/
EstimateParametersHMRFEM.m
42 lines (37 loc) · 1.46 KB
/
EstimateParametersHMRFEM.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
%% Îöåíèâàåò ïàðàìåòðû beta, mu, kappa ïî äàííûì è ìàòðèöå àïîñòåðèîðíûõ âåðîÿòíîñòåé
%---input---------------------------------------------------------
% X: 4-õ ìåðíàÿ ìàòðèöà èñõîäíûõ äàííûõ
% posterior: àïîñòåðèîðíûå âåðîÿòíîñòè, LxN
% k: êîëè÷åñòâî êëàññîâ ñåãìåíòàöèè
% p: ðàçìåðíîñòü äàííûõ
% beta: òåêóùåå çíà÷åíèå ïàðàìåòðà ìîäåëè Ïîòòñà
% mu: òåêóùàÿ ìàòðèöà ïàðàìåòðîâ äëÿ vMF
% kappa: òåêóùèé âåêòîð ïàðàìåòðîâ äëÿ vMF
%---output--------------------------------------------------------
% beta: íîâîå çíà÷åíèå ïàðàìåòðà, ñêàëÿð
% mu: íîâûå çíà÷åíèÿ ïàðàìåòðîâ, ìàòðèöà LxP
% kappa: íîâûå çíà÷åíèÿ ïàðàìåòðîâ, âåêòîð 1xL
function [beta, mu, kappa] = EstimateParametersHMRFEM(X, Y, posterior, k, p, beta, mu, kappa)
for l=1:k
R = posterior(l,Y==l) * X(Y==l,:);
Rlen = sqrt(sum(R .^ 2));
if Rlen ~= 0
mu(l, :) = R / Rlen;
end
eqn = @(x) besseli(p/2, x, 1)/besseli(p/2-1,x, 1) * sum(posterior(l,Y==l)) - Rlen;
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