diff --git a/nonlinear-filters/src/conditional_filter_proposal.m b/nonlinear-filters/src/conditional_filter_proposal.m index 5137b83f7..fe9278c48 100644 --- a/nonlinear-filters/src/conditional_filter_proposal.m +++ b/nonlinear-filters/src/conditional_filter_proposal.m @@ -1,4 +1,4 @@ -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 % for the state errors and the Kalman filter @@ -42,6 +42,7 @@ persistent init_flag2 mf0 mf1 persistent number_of_state_variables number_of_observed_variables persistent number_of_structural_innovations +flag=0 ; % Set local state space model (first-order approximation). ghx = ReducedForm.ghx; ghu = ReducedForm.ghu; @@ -118,15 +119,31 @@ else Error = obs - PredictedObservedMean ; StateVectorMean = PredictedStateMean + KalmanFilterGain*Error ; 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 Weights = SampleWeights.*probability2(zeros(number_of_observed_variables,1),chol(PredictedObservedVariance)',Error) ; end end -PredictedStateVarianceSquareRoot = chol(PredictedStateVariance + eye(number_of_state_variables)*1e-6)' ; ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot,2),1)+StateVectorMean ; 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) ; Posterior = probability2(StateVectorMean,StateVectorVarianceSquareRoot,ProposalStateVector) ; Likelihood = probability2(obs,H_lower_triangular_cholesky,measurement_equations(ProposalStateVector,ReducedForm,ThreadsOptions)) ; diff --git a/nonlinear-filters/src/conditional_particle_filter.m b/nonlinear-filters/src/conditional_particle_filter.m index 6216cf15c..9ba069a19 100644 --- a/nonlinear-filters/src/conditional_particle_filter.m +++ b/nonlinear-filters/src/conditional_particle_filter.m @@ -103,8 +103,13 @@ StateParticles = bsxfun(@plus,StateVectorVarianceSquareRoot*randn(state_variance SampleWeights = ones(1,number_of_particles)/number_of_particles ; for t=1:sample_size 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) ; + if flag==1 + LIK=-Inf; + lik(t)=-Inf; + return + end end SumSampleWeights = sum(SampleWeights) ; lik(t) = log(SumSampleWeights) ; @@ -116,4 +121,4 @@ for t=1:sample_size end end -LIK = -sum(lik(start:end)); \ No newline at end of file +LIK = -sum(lik(start:end));