diff --git a/matlab/nonlinear-filters/nonlinear_kalman_filter.m b/matlab/nonlinear-filters/nonlinear_kalman_filter.m
index 354100c2f04be2f90135e25f95e2546a22a828fb..551345c1acd2aade60c6658b08b5143994a9ba58 100644
--- a/matlab/nonlinear-filters/nonlinear_kalman_filter.m
+++ b/matlab/nonlinear-filters/nonlinear_kalman_filter.m
@@ -123,6 +123,8 @@ StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
 lik  = NaN(sample_size,1);
 LIK  = NaN;
 
+n = rows(Y);
+
 for t=1:sample_size
     xbar = [StateVectorMean ; zeros(number_of_structural_innovations,1) ] ;
     sqr_Px = [StateVectorVarianceSquareRoot zeros(number_of_state_variables,number_of_structural_innovations);
@@ -186,7 +188,10 @@ for t=1:sample_size
             return
         end
     end
-    lik(t) = log( sum(probability2(Y(:,t),H_lower_triangular_cholesky,tmp(mf1,:)).*weights,1) ) ;
+    PredictedObservedVarianceSquareRoot = chol(PredictedObservedVariance);
+    PredictedObservedInverseVarianceSquareRoot = PredictedObservedVarianceSquareRoot\eye(n); % Inverse of the Cholesky -> inv(PredictedObservedVariance) = A*A'
+    StandardPredictionError = PredictedObservedInverseVarianceSquareRoot'*PredictionError;
+    lik(t) = -.5*n*log(2*pi) + sum(log(diag(PredictedObservedInverseVarianceSquareRoot))) - .5*dot(StandardPredictionError, StandardPredictionError);
 end
 
 LIK = -sum(lik(start:end));
diff --git a/meson.build b/meson.build
index 89abf1f985eb5c79078a49a47d65a3dcd9ddb117..1f162f7627dd87dd3632c391511936d9a0813fe1 100644
--- a/meson.build
+++ b/meson.build
@@ -1770,6 +1770,9 @@ mod_and_m_tests = [
   { 'test' : [ 'particle/noisy_ar1.mod',
                'particle/noisy_ar1_nlf_sis.mod'],
     'extra' : [ 'particle/noisy_ar1_data.m'] },
+  { 'test' : [ 'particle/noisy_ar1.mod',
+               'particle/noisy_ar1_nlf_nlkf.mod'],
+    'extra' : [ 'particle/noisy_ar1_data.m'] },
   { 'test' : [ 'particle/dsge_base2.mod' ],
     'extra' : [ 'particle/risky.m',
                 'particle/extreme.m',