computeGammaHMM.m 2.89 KB
Newer Older
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 80
function [GAMMA, ALPHA] = computeGammaHMM(s, model)
% Compute Gamma probability in HMM by using forward and backward variables.
%
% Writing code takes time. Polishing it and making it available to others takes longer! 
% If some parts of the code were useful for your research of for a better understanding 
% of the algorithms, please reward the authors by citing the related publications, 
% and consider making your own research available in this way.
%
% @article{Rozo16Frontiers,
%   author="Rozo, L. and Silv\'erio, J. and Calinon, S. and Caldwell, D. G.",
%   title="Learning Controllers for Reactive and Proactive Behaviors in Human-Robot Collaboration",
%   journal="Frontiers in Robotics and {AI}",
%   year="2016",
%   month="June",
%   volume="3",
%   number="30",
%   pages="1--11",
%   doi="10.3389/frobt.2016.00030"
% }
% 
% Copyright (c) 2015 Idiap Research Institute, http://idiap.ch/
% Written by Sylvain Calinon, http://calinon.ch/
% 
% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
% 
% PbDlib is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 3 as
% published by the Free Software Foundation.
% 
% PbDlib is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with PbDlib. If not, see <http://www.gnu.org/licenses/>.


%Initialization of the parameters
nbSamples = length(s);
Data = [];
for n=1:nbSamples
	Data = [Data s(n).Data];
	s(n).nbData = size(s(n).Data,2);
end

for n=1:nbSamples
	%Emission probabilities
	for i=1:model.nbStates
		s(n).B(i,:) = model.Priors(i) * gaussPDF(s(n).Data, model.Mu(:,i), model.Sigma(:,:,i));
	end
	%Forward variable ALPHA
	s(n).ALPHA(:,1) = model.StatesPriors .* s(n).B(:,1);
	%Scaling to avoid underflow issues
	s(n).c(1) = 1 / sum(s(n).ALPHA(:,1)+realmin);
	s(n).ALPHA(:,1) = s(n).ALPHA(:,1) * s(n).c(1);
	for t=2:s(n).nbData
		s(n).ALPHA(:,t) = (s(n).ALPHA(:,t-1)'*model.Trans)' .* s(n).B(:,t); 
		%Scaling to avoid underflow issues
		s(n).c(t) = 1 / sum(s(n).ALPHA(:,t)+realmin);
		s(n).ALPHA(:,t) = s(n).ALPHA(:,t) * s(n).c(t);
	end
	%Backward variable BETA
	s(n).BETA(:,s(n).nbData) = ones(model.nbStates,1) * s(n).c(end); %Rescaling
	for t=s(n).nbData-1:-1:1
		s(n).BETA(:,t) = model.Trans * (s(n).BETA(:,t+1) .* s(n).B(:,t+1));
		s(n).BETA(:,t) = min(s(n).BETA(:,t) * s(n).c(t), realmax); %Rescaling
	end
	%Intermediate variable GAMMA
	s(n).GAMMA = (s(n).ALPHA.*s(n).BETA) ./ repmat(sum(s(n).ALPHA.*s(n).BETA)+realmin, model.nbStates, 1); 
end

%Concatenation of GAMMAs
GAMMA=[]; ALPHA=[]; 
for n=1:nbSamples
	GAMMA = [GAMMA s(n).GAMMA];
	ALPHA = [ALPHA s(n).ALPHA];
end