From cd350a6f578a6e343965533740ed192e3f9920c7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Argos=29?=
<stepan@adjemian.eu>
Date: Mon, 14 Apr 2025 20:23:31 +0200
Subject: [PATCH] [WIP] Fix auxiliary particle filter.
---
.../nonlinear-filters/auxiliary_particle_filter.m | 13 ++++++++-----
1 file changed, 8 insertions(+), 5 deletions(-)
diff --git a/matlab/nonlinear-filters/auxiliary_particle_filter.m b/matlab/nonlinear-filters/auxiliary_particle_filter.m
index 74354dc69e..803cff035c 100644
--- a/matlab/nonlinear-filters/auxiliary_particle_filter.m
+++ b/matlab/nonlinear-filters/auxiliary_particle_filter.m
@@ -76,13 +76,14 @@ H = ReducedForm.H;
StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2);
-Q_lower_triangular_cholesky = chol(Q)';
+Q_lower_triangular_cholesky = chol(Q, 'lower');
+H_lower_triangular_cholesky = chol(H, 'lower');
% Set seed for randn().
options_=set_dynare_seed_local_options(options_,'default');
% Initialization of the likelihood.
-const_lik = log(2*pi)*number_of_observed_variables+log(det(H));
+const_lik = log(2*pi)*number_of_observed_variables + sum(log(diag(H_lower_triangular_cholesky)));
lik = NaN(sample_size,1);
LIK = NaN;
@@ -140,13 +141,15 @@ for t=1:sample_size
end
end
end
- PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
- z = sum(PredictionError.*(H\PredictionError),1) ;
+ StandardPredictionError = H_lower_triangular_cholesky\bsxfun(@minus,Y(:,t), tmp(mf1,:));
+ lambda = weights.*exp(-.5*(const_lik + sum(StandardPredictionError.*StandardPredictionError, 1)));
+ LAMBDA = lambda/sum(lambda);
+ indx = resample(0, LAMBDA, ParticleOptions);
% tau_tilde = weights.*(tpdf(z,3*ones(size(z)))+1e-99) ;
ddl = 3 ;
tau_tilde = weights.*(exp(gammaln((ddl + 1) / 2) - gammaln(ddl/2))./(sqrt(ddl*pi).*(1 + (z.^2)./ddl).^((ddl + 1)/2))+1e-99) ;
tau_tilde = tau_tilde/sum(tau_tilde) ;
- indx = resample(0,tau_tilde',ParticleOptions);
+
if pruning
yhat_ = yhat_(:,indx) ;
end
--
GitLab