Commit aad6c297 authored by Frédéric Karamé's avatar Frédéric Karamé
Browse files

Reintroduction of pruning and of options to build the sufficient statistics in...

Reintroduction of pruning and of options to build the sufficient statistics in the auxiliary particle filter
parent 473b2b8b
...@@ -98,21 +98,32 @@ else ...@@ -98,21 +98,32 @@ else
end end
StateVectors = tmp(mf0,:) ; StateVectors = tmp(mf0,:) ;
%[nodes,nodes_weights] = spherical_radial_sigma_points(number_of_structural_innovations) ; if ParticleOptions.proposal_approximation.cubature
[nodes,nodes_weights,nodes_weights_c] = unscented_sigma_points(number_of_structural_innovations,ParticleOptions); [nodes,nodes_weights] = spherical_radial_sigma_points(number_of_structural_innovations) ;
nodes_weights = ones(size(nodes,1),1)*nodes_weights ;
elseif ParticleOptions.proposal_approximation.unscented
[nodes,nodes_weights,nodes_weights_c] = unscented_sigma_points(number_of_structural_innovations,ParticleOptions);
else
error('Estimation: This approximation for the proposal is not implemented or unknown!')
end
nodes = Q_lower_triangular_cholesky*nodes ; nodes = Q_lower_triangular_cholesky*nodes ;
for t=1:sample_size for t=1:sample_size
yhat = bsxfun(@minus,StateVectors,state_variables_steady_state); yhat = bsxfun(@minus,StateVectors,state_variables_steady_state);
% if pruning if pruning
% yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state); yhat_ = bsxfun(@minus,StateVectors_,state_variables_steady_state);
% [tmp, tmp_] = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2); tmp = 0 ;
% else tmp_ = 0 ;
% tmp = local_state_space_iteration_2(yhat,zeros(number_of_structural_innovations,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); for i=1:size(nodes)
% end [tmp1, tmp1_] = local_state_space_iteration_2(yhat,nodes(i,:)*ones(1,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,yhat_,steadystate,ThreadsOptions.local_state_space_iteration_2);
tmp = 0 ; tmp = tmp + nodes_weights(i)*tmp1 ;
for i=1:size(nodes) tmp_ = tmp_ + nodes_weights(i)*tmp1_ ;
tmp = tmp + nodes_weights*local_state_space_iteration_2(yhat,nodes(i,:)*ones(1,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); end
else
tmp = 0 ;
for i=1:size(nodes)
tmp = tmp + nodes_weights(i)*local_state_space_iteration_2(yhat,nodes(i,:)*ones(1,number_of_particles),ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
end
end end
PredictedObservedMean = weights*(tmp(mf1,:)'); PredictedObservedMean = weights*(tmp(mf1,:)');
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
...@@ -128,7 +139,6 @@ for t=1:sample_size ...@@ -128,7 +139,6 @@ for t=1:sample_size
yhat_ = yhat_(:,indx) ; yhat_ = yhat_(:,indx) ;
end end
yhat = yhat(:,indx) ; yhat = yhat(:,indx) ;
%wtilde = wtilde(indx) ;
factor = weights(indx)./tau_tilde(indx) ; factor = weights(indx)./tau_tilde(indx) ;
epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles); epsilon = Q_lower_triangular_cholesky*randn(number_of_structural_innovations,number_of_particles);
if pruning if pruning
...@@ -143,7 +153,6 @@ for t=1:sample_size ...@@ -143,7 +153,6 @@ for t=1:sample_size
dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean); dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
PredictedObservedVariance = (dPredictedObservedMean*dPredictedObservedMean')/number_of_particles + H; PredictedObservedVariance = (dPredictedObservedMean*dPredictedObservedMean')/number_of_particles + H;
lnw = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))); lnw = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1)));
%wtilde = lnw./wtilde;
wtilde = lnw.*factor ; wtilde = lnw.*factor ;
weights = wtilde/sum(wtilde); weights = wtilde/sum(wtilde);
end end
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment