Commit b8db5a23 authored by fred's avatar fred
Browse files

latest version ; still buggy

git-svn-id: https://www.dynare.org/svn/dynare/trunk@3232 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 9db0e3a8
......@@ -45,7 +45,7 @@ persistent sample_size number_of_state_variables number_of_observed_variables nu
% Set defaults.
if (nargin<4) || (nargin==4 && isempty(number_of_particles))
number_of_particles = 1000 ;
number_of_particles = 10 ;
end
if nargin==2 || isempty(start)
start = 1;
......@@ -96,6 +96,7 @@ randn('state',seed);
const_lik = log(2*pi)*number_of_observed_variables;
lik = NaN(sample_size,1);
nb_obs_resamp = 0 ;
w = ones(number_of_particles,1) ;
for t=1:sample_size
PredictedState = zeros(number_of_particles,number_of_state_variables);
PredictionError = zeros(number_of_particles,number_of_observed_variables);
......@@ -133,48 +134,47 @@ for t=1:sample_size
end
PredictedObservedVariance = PredictedObservedVariance + H;
iPredictedObservedVariance = inv(PredictedObservedVariance);
lnw = - 0.5*(const_lik + log(det(PredictedObservedVariance)) + sum(PredictionError'*iPredictedObservedVariance.*PredictionError',2)) ;
lnw = -0.5*(const_lik + log(det(PredictedObservedVariance)) + sum((PredictionError*iPredictedObservedVariance).*PredictionError,2)) ;
%bidouille numrique Schorfheide
dfac = max(lnw);
wtilde = w.*exp(lnw - dfac) ;
% vraisemblance de l'observation
lik(t) = log(mean(wtilde)) + dfac ;
clear (PredictionError) ;
clear (lnw) ;
%clear (PredictionError) ;
%clear (lnw) ;
% calcul des poids
w = wtilde/sum(wtilde) ;
clear (wtilde) ;
%clear (wtilde) ;
%update
Neff = 1/sum(w^2) ;
if Neff>number_of_particules %no resampling
Neff = 1/sum(w.*w) ;
if Neff>number_of_particles %no resampling
StateUpdated = PredictedState ;
clear (PredictedState) ;
%clear (PredictedState) ;
w = number_of_particles*w ;
else %resampling
nb_obs_resamp = nb_obs_resamp+1 ;
%kill the smallest particles before resampling :! facultatif ?
to_kill = [w PredictedState] ;
to_kill = delif(to_kill,w<(1/number_of_particules)*1E-12);%%
[n,m] = size(to_kill) ;
w = to_kill(:,1) ;
PredictedState = to_kill(:,2:m) ;
clear (to_kill) ;
if number_of_particles neq n
'Elimination de '; number_of_particles - n ; ' particules l''observation ';t ;
end
%to_kill = [w PredictedState] ;
%to_kill = delif(to_kill,w<(1/number_of_particules)*1E-12);%%
%[n,m] = size(to_kill) ;
%w = to_kill(:,1) ;
%PredictedState = to_kill(:,2:m) ;
%clear (to_kill) ;
%if number_of_particles neq n
% 'Elimination de '; number_of_particles - n ; ' particules l''observation ';t ;
%end
%fin de kill
%remise l'chelle des poids sur les particules restantes
w = cumsum( w/sum(w) );
%w = cumsum( w/sum(w) );
%Rchantillonage systmatique
rnduvec = ( (1:number_of_particles)-1+rand )/number_of_particles ;
selind = (n - sum( w > rnduvec' ) + 1)'; % problme de mmoire car w .> rnduvec' trs grande !
clear (rnduvec) ;
StateUpdated = PredictedState(selind,:) ;
clear (selind) ;
%rnduvec = ( (1:number_of_particles)-1+rand )/number_of_particles ;
%selind = (number_of_particles - sum( w > rnduvec ) + 1)'; % problme de mmoire car w .> rnduvec' trs grande !
%clear (rnduvec) ;
%StateUpdated = PredictedState(selind,:) ;
%clear (selind) ;
% initialize
selind = zeros(1,number_of_particles);
selind = zeros(number_of_particles,1);
% construct CDF
c = cumsum(w);
% draw a starting point
......@@ -190,7 +190,7 @@ for t=1:sample_size
selind(i) = j;
end
StateUpdated = PredictedState(selind,:);
w = ones(number_of_particules,1) ;
w = ones(number_of_particles,1) ;
end
end
LIK = -sum(lik(start:end));
\ No newline at end of file
Supports Markdown
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