diff --git a/m_fcts/EM_GMM.m b/m_fcts/EM_GMM.m
index 8e8339613c1742253c7740ff735647692e6d48fe..710d182b2fbbb041e60973617ae951189c6a0dcb 100644
--- a/m_fcts/EM_GMM.m
+++ b/m_fcts/EM_GMM.m
@@ -59,7 +59,7 @@ function [L, GAMMA] = computeGamma(Data, model)
   for i=1:model.nbStates
 		L(i,:) = model.Priors(i) * gaussPDF(Data, model.Mu(:,i), model.Sigma(:,:,i));
 	end
-	GAMMA = L ./ repmat(sum(L,1),model.nbStates,1);
+	GAMMA = L ./ repmat(sum(L,1)+realmin,model.nbStates,1);
 end
 
 
diff --git a/m_fcts/GMR.m b/m_fcts/GMR.m
index 5e6e4ccce5646a5afdb33d6bf2699d7d04fd1926..956de0355d430f6370d4b3cdf65ad893c9409fa8 100644
--- a/m_fcts/GMR.m
+++ b/m_fcts/GMR.m
@@ -20,18 +20,18 @@ function [expData, expSigma, H] = GMR(model, DataIn, in, out)
 nbData = size(DataIn,2);
 nbVarOut = length(out);
 
+diagRegularizationFactor = 1E-8;
+
 MuTmp = zeros(nbVarOut,model.nbStates);
 expData = zeros(nbVarOut,nbData);
 expSigma = zeros(nbVarOut,nbVarOut,nbData);
-expSigma2 = zeros(nbVarOut,nbVarOut,nbData);
-expSigma3 = zeros(nbVarOut,nbVarOut,nbData);
 for t=1:nbData
 	%Compute activation weight
 	%See Eq. (3.0.5) in doc/TechnicalReport.pdf
 	for i=1:model.nbStates
 		H(i,t) = model.Priors(i) * gaussPDF(DataIn(:,t), model.Mu(in,i), model.Sigma(in,in,i));
 	end
-	H(:,t) = H(:,t)/sum(H(:,t));
+	H(:,t) = H(:,t)/sum(H(:,t)+realmin);
 	%Compute expected conditional means
 	%See Eq. (3.0.8) in doc/TechnicalReport.pdf
 	for i=1:model.nbStates
@@ -44,6 +44,6 @@ for t=1:nbData
 		SigmaTmp = model.Sigma(out,out,i) - model.Sigma(out,in,i)/model.Sigma(in,in,i) * model.Sigma(in,out,i);
 		expSigma(:,:,t) = expSigma(:,:,t) + H(i,t) * (SigmaTmp + MuTmp(:,i)*MuTmp(:,i)');
 	end
-	expSigma(:,:,t) = expSigma(:,:,t) - expData(:,t)*expData(:,t)'; 
+	expSigma(:,:,t) = expSigma(:,:,t) - expData(:,t)*expData(:,t)' + eye(nbVarOut) * diagRegularizationFactor; 
 end