stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fix a bug in likelihood calculation.

rm-particles^2
Frédéric Karamé 2015-10-02 16:05:26 +02:00
parent cdc7f6ddf9
commit 04ad104bfb
6 changed files with 30 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
src/auxiliary_particle_filter.m
View File

@ -64,7 +64,7 @@ end
% Get initial condition for the state vector. % Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean; StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)'; StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';%reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2); state_variance_rank = size(StateVectorVarianceSquareRoot,2);
Q_lower_triangular_cholesky = chol(Q)'; Q_lower_triangular_cholesky = chol(Q)';
if pruning if pruning
@ -76,7 +76,7 @@ end
set_dynare_seed('default'); set_dynare_seed('default');
% Initialization of the likelihood. % Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables; const_lik = log(2*pi)*number_of_observed_variables +log(det(H));
lik = NaN(sample_size,1); lik = NaN(sample_size,1);
LIK = NaN; LIK = NaN;
@ -125,11 +125,11 @@ for t=1:sample_size
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); 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 end
PredictedObservedMean = weights*(tmp(mf1,:)'); %PredictedObservedMean = weights*(tmp(mf1,:)');
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean'); %dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean');
PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' +H; %PredictedObservedVariance = bsxfun(@times,weights,dPredictedObservedMean)*dPredictedObservedMean' +H;
wtilde = exp(-.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1))) ; wtilde = exp(-.5*(const_lik+sum(PredictionError.*(H\PredictionError),1))) ;
tau_tilde = weights.*wtilde ; tau_tilde = weights.*wtilde ;
sum_tau_tilde = sum(tau_tilde) ; sum_tau_tilde = sum(tau_tilde) ;
lik(t) = log(sum_tau_tilde) ; lik(t) = log(sum_tau_tilde) ;
@ -148,11 +148,11 @@ for t=1:sample_size
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
end end
StateVectors = tmp(mf0,:); StateVectors = tmp(mf0,:);
PredictedObservedMean = mean(tmp(mf1,:),2); %PredictedObservedMean = mean(tmp(mf1,:),2);
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
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+sum(PredictionError.*(H\PredictionError),1)));
wtilde = lnw.*factor ; wtilde = lnw.*factor ;
weights = wtilde/sum(wtilde); weights = wtilde/sum(wtilde);
end end

2
src/conditional_filter_proposal.m
View File

@ -113,7 +113,7 @@ else
StateVectorMean = PredictedStateMean + KalmanFilterGain*(obs - PredictedObservedMean); StateVectorMean = PredictedStateMean + KalmanFilterGain*(obs - PredictedObservedMean);
StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain'; StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance'); StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance');
StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)'; StateVectorVarianceSquareRoot = chol(StateVectorVariance)';%reduced_rank_cholesky(StateVectorVariance)';
end end
ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot,2),1)+StateVectorMean ; ProposalStateVector = StateVectorVarianceSquareRoot*randn(size(StateVectorVarianceSquareRoot,2),1)+StateVectorMean ;

12
src/conditional_particle_filter.m
View File

@ -57,9 +57,9 @@ function [LIK,lik] = conditional_particle_filter(ReducedForm,Y,start,ParticleOpt
% AUTHOR(S) frederic DOT karame AT univ DASH lemans DOT fr % AUTHOR(S) frederic DOT karame AT univ DASH lemans DOT fr
% stephane DOT adjemian AT univ DASH lemans DOT fr % stephane DOT adjemian AT univ DASH lemans DOT fr
persistent init_flag mf0 mf1 persistent init_flag mf1
persistent number_of_particles persistent number_of_particles
persistent sample_size number_of_state_variables number_of_observed_variables persistent sample_size number_of_observed_variables
% Set default % Set default
if isempty(start) if isempty(start)
@ -68,10 +68,10 @@ end
% Set persistent variables. % Set persistent variables.
if isempty(init_flag) if isempty(init_flag)
mf0 = ReducedForm.mf0; %mf0 = ReducedForm.mf0;
mf1 = ReducedForm.mf1; mf1 = ReducedForm.mf1;
sample_size = size(Y,2); sample_size = size(Y,2);
number_of_state_variables = length(mf0); %number_of_state_variables = length(mf0);
number_of_observed_variables = length(mf1); number_of_observed_variables = length(mf1);
init_flag = 1; init_flag = 1;
number_of_particles = ParticleOptions.number_of_particles ; number_of_particles = ParticleOptions.number_of_particles ;
@ -84,14 +84,14 @@ if isempty(H)
H = 0; H = 0;
H_lower_triangular_cholesky = 0; H_lower_triangular_cholesky = 0;
else else
H_lower_triangular_cholesky = reduced_rank_cholesky(H)'; H_lower_triangular_cholesky = chol(H)'; %reduced_rank_cholesky(H)';
end end
% Get initial condition for the state vector. % Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean; StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)'; StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
state_variance_rank = size(StateVectorVarianceSquareRoot,2); state_variance_rank = size(StateVectorVarianceSquareRoot,2);
Q_lower_triangular_cholesky = reduced_rank_cholesky(Q)'; Q_lower_triangular_cholesky = chol(Q)'; %reduced_rank_cholesky(Q)';
% Set seed for randn(). % Set seed for randn().
set_dynare_seed('default'); set_dynare_seed('default');

7
src/gaussian_densities.m
View File

@ -43,9 +43,10 @@ prior = probability2(st_t_1,sqr_Pss_t_t_1,particles) ;
% likelihood % likelihood
yt_t_1_i = measurement_equations(particles,ReducedForm,ThreadsOptions) ; yt_t_1_i = measurement_equations(particles,ReducedForm,ThreadsOptions) ;
eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ; eta_t_i = bsxfun(@minus,obs,yt_t_1_i)' ;
yt_t_1 = sum(yt_t_1_i*weigths1,2) ; %yt_t_1 = sum(yt_t_1_i*weigths1,2) ;
tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ; %tmp = bsxfun(@minus,yt_t_1_i,yt_t_1) ;
Pyy = bsxfun(@times,weigths2',tmp)*tmp' + H ; %Pyy = bsxfun(@times,weigths2',tmp)*tmp' + H ;
Pyy = H ;
sqr_det = sqrt(det(Pyy)) ; sqr_det = sqrt(det(Pyy)) ;
foo = (eta_t_i/Pyy).*eta_t_i ; foo = (eta_t_i/Pyy).*eta_t_i ;
likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ; likelihood = exp(-0.5*sum(foo,2))/(normconst*sqr_det) + 1e-99 ;

2
src/gaussian_filter_bank.m
View File

@ -118,5 +118,5 @@ else
StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError; StateVectorMean = PredictedStateMean + KalmanFilterGain*PredictionError;
StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain'; StateVectorVariance = PredictedStateVariance - KalmanFilterGain*PredictedObservedVariance*KalmanFilterGain';
StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance'); StateVectorVariance = .5*(StateVectorVariance+StateVectorVariance');
StateVectorVarianceSquareRoot = reduced_rank_cholesky(StateVectorVariance)'; StateVectorVarianceSquareRoot = chol(StateVectorVariance)'; %reduced_rank_cholesky(StateVectorVariance)';
end end

15
src/sequential_importance_particle_filter.m
View File

@ -66,12 +66,12 @@ if isempty(H)
end end
% Initialization of the likelihood. % Initialization of the likelihood.
const_lik = log(2*pi)*number_of_observed_variables; const_lik = log(2*pi)*number_of_observed_variables +log(det(H)) ;
lik = NaN(sample_size,1); lik = NaN(sample_size,1);
% Get initial condition for the state vector. % Get initial condition for the state vector.
StateVectorMean = ReducedForm.StateVectorMean; StateVectorMean = ReducedForm.StateVectorMean;
StateVectorVarianceSquareRoot = reduced_rank_cholesky(ReducedForm.StateVectorVariance)'; StateVectorVarianceSquareRoot = chol(ReducedForm.StateVectorVariance)';%reduced_rank_cholesky(ReducedForm.StateVectorVariance)';
if pruning if pruning
StateVectorMean_ = StateVectorMean; StateVectorMean_ = StateVectorMean;
StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot; StateVectorVarianceSquareRoot_ = StateVectorVarianceSquareRoot;
@ -103,12 +103,13 @@ for t=1:sample_size
else else
tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2); tmp = local_state_space_iteration_2(yhat,epsilon,ghx,ghu,constant,ghxx,ghuu,ghxu,ThreadsOptions.local_state_space_iteration_2);
end end
PredictedObservedMean = tmp(mf1,:)*transpose(weights); %PredictedObservedMean = tmp(mf1,:)*transpose(weights);
PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:)); PredictionError = bsxfun(@minus,Y(:,t),tmp(mf1,:));
dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean); %dPredictedObservedMean = bsxfun(@minus,tmp(mf1,:),PredictedObservedMean);
PredictedObservedVariance = bsxfun(@times,dPredictedObservedMean,weights)*dPredictedObservedMean' + H; %PredictedObservedVariance = bsxfun(@times,dPredictedObservedMean,weights)*dPredictedObservedMean' + H;
if rcond(PredictedObservedVariance) > 1e-16 %PredictedObservedVariance = H;
lnw = -.5*(const_lik+log(det(PredictedObservedVariance))+sum(PredictionError.*(PredictedObservedVariance\PredictionError),1)); if rcond(H) > 1e-16
lnw = -.5*(const_lik+sum(PredictionError.*(H\PredictionError),1));
else else
LIK = NaN; LIK = NaN;
return return