diff --git a/src/gaussian_mixture_densities.m b/src/gaussian_mixture_densities.m
index c826059345220f89b85bdd0908c17d9daba76c4d..cd4f6b58e61326e3fdc2d34b9fa161acee370952 100644
--- a/src/gaussian_mixture_densities.m
+++ b/src/gaussian_mixture_densities.m
@@ -46,11 +46,12 @@ prior = prior' ;
 proposal = proposal' ;
 % Compute the density of the current observation conditionally to each particle
 yt_t_1_i = measurement_equations(StateParticles,ReducedForm,ThreadsOptions) ;
-eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ;
-yt_t_1 = sum(yt_t_1_i*weigths1,2) ;
-tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;
-Pyy = bsxfun(@times,weigths2',tmp)*tmp' + H ;
-sqr_det = sqrt(det(Pyy)) ;
-foo = (eta_t_i/Pyy).*eta_t_i ;
-likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ;
+%eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ;
+%yt_t_1 = sum(yt_t_1_i*weigths1,2) ;
+%tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;
+%Pyy = bsxfun(@times,weigths2',tmp)*tmp' + H ;
+%sqr_det = sqrt(det(Pyy)) ;
+%foo = (eta_t_i/Pyy).*eta_t_i ;
+%likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ;
+likelihood = probability2(obs,sqrt(H),yt_t_1_i) ;
 IncrementalWeights = likelihood.*prior./proposal ;