Commit 43615ce4 authored by Frédéric Karamé's avatar Frédéric Karamé

Add flags to deal with errors on Cholesky matrices in CPF filter.

parent 1f08164b
function [ProposalStateVector,Weights] = conditional_filter_proposal(ReducedForm,obs,StateVectors,SampleWeights,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions,normconst2) function [ProposalStateVector,Weights,flag] = conditional_filter_proposal(ReducedForm,obs,StateVectors,SampleWeights,Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions,normconst2)
% %
% Computes the proposal for each past particle using Gaussian approximations % Computes the proposal for each past particle using Gaussian approximations
% for the state errors and the Kalman filter % for the state errors and the Kalman filter
...@@ -42,6 +42,7 @@ persistent init_flag2 mf0 mf1 ...@@ -42,6 +42,7 @@ persistent init_flag2 mf0 mf1
persistent number_of_state_variables number_of_observed_variables persistent number_of_state_variables number_of_observed_variables
persistent number_of_structural_innovations persistent number_of_structural_innovations
flag=0 ;
% Set local state space model (first-order approximation). % Set local state space model (first-order approximation).
ghx = ReducedForm.ghx; ghx = ReducedForm.ghx;
ghu = ReducedForm.ghu; ghu = ReducedForm.ghu;
...@@ -118,15 +119,31 @@ else ...@@ -118,15 +119,31 @@ else
Error = obs - PredictedObservedMean ; Error = obs - PredictedObservedMean ;
StateVectorMean = PredictedStateMean + KalmanFilterGain*Error ; StateVectorMean = PredictedStateMean + KalmanFilterGain*Error ;
StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain'; StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
StateVectorVarianceSquareRoot = chol(StateVectorVariance + eye(number_of_state_variables)*1e-6)' ; StateVectorVariance = 0.5*(StateVectorVariance+StateVectorVariance');
%StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)';%chol(StateVectorVariance + eye(number_of_state_variables)*1e-6)' ;
[StateVectorVarianceSquareRoot,p] = chol(StateVectorVariance,'lower') ;
if p
flag=1;
ProposalStateVector = zeros(number_of_state_variables,1) ;
Weights = 0.0 ;
return
end
if ParticleOptions.cpf_weights_method.amisanotristani if ParticleOptions.cpf_weights_method.amisanotristani
Weights = SampleWeights.*probability2(zeros(number_of_observed_variables,1),chol(PredictedObservedVariance)',Error) ; Weights = SampleWeights.*probability2(zeros(number_of_observed_variables,1),chol(PredictedObservedVariance)',Error) ;
end end
end end
PredictedStateVarianceSquareRoot = chol(PredictedStateVariance + eye(number_of_state_variables)*1e-6)' ;
ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot,2),1)+StateVectorMean ; ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot,2),1)+StateVectorMean ;
if ParticleOptions.cpf_weights_method.murrayjonesparslow if ParticleOptions.cpf_weights_method.murrayjonesparslow
PredictedStateVariance = 0.5*(PredictedStateVariance+PredictedStateVariance');
%PredictedStateVarianceSquareRoot = reduced_rank_cholesky(PredictedStateVariance)';%chol(PredictedStateVariance + eye(number_of_state_variables)*1e-6)' ;
[PredictedStateVarianceSquareRoot,p] = chol(PredictedStateVariance,'lower') ;
if p
flag=1;
ProposalStateVector = zeros(number_of_state_variables,1) ;
Weights = 0.0 ;
return
end
Prior = probability2(PredictedStateMean,PredictedStateVarianceSquareRoot,ProposalStateVector) ; Prior = probability2(PredictedStateMean,PredictedStateVarianceSquareRoot,ProposalStateVector) ;
Posterior = probability2(StateVectorMean,StateVectorVarianceSquareRoot,ProposalStateVector) ; Posterior = probability2(StateVectorMean,StateVectorVarianceSquareRoot,ProposalStateVector) ;
Likelihood = probability2(obs,H_lower_triangular_cholesky,measurement_equations(ProposalStateVector,ReducedForm,ThreadsOptions)) ; Likelihood = probability2(obs,H_lower_triangular_cholesky,measurement_equations(ProposalStateVector,ReducedForm,ThreadsOptions)) ;
......
...@@ -103,8 +103,13 @@ StateParticles = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance ...@@ -103,8 +103,13 @@ StateParticles = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance
SampleWeights = ones(1,number_of_particles)/number_of_particles ; SampleWeights = ones(1,number_of_particles)/number_of_particles ;
for t=1:sample_size for t=1:sample_size
for i=1:number_of_particles for i=1:number_of_particles
[StateParticles(:,i),SampleWeights(i)] = ... [StateParticles(:,i),SampleWeights(i),flag] = ...
conditional_filter_proposal(ReducedForm,Y(:,t),StateParticles(:,i),SampleWeights(i),Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions,normconst2) ; conditional_filter_proposal(ReducedForm,Y(:,t),StateParticles(:,i),SampleWeights(i),Q_lower_triangular_cholesky,H_lower_triangular_cholesky,H,ParticleOptions,ThreadsOptions,normconst2) ;
if flag==1
LIK=-Inf;
lik(t)=-Inf;
return
end
end end
SumSampleWeights = sum(SampleWeights) ; SumSampleWeights = sum(SampleWeights) ;
lik(t) = log(SumSampleWeights) ; lik(t) = log(SumSampleWeights) ;
...@@ -116,4 +121,4 @@ for t=1:sample_size ...@@ -116,4 +121,4 @@ for t=1:sample_size
end end
end end
LIK = -sum(lik(start:end)); LIK = -sum(lik(start:end));
\ No newline at end of file
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