Commit a73e4d94 by Frédéric Karamé

### Add the possibility of Gaussian-Mixture Particle Filter without resampling.

parent a6a3e1b1
 function [StateMu,StateSqrtP,StateWeights] = fit_gaussian_mixture(X,StateMu,StateSqrtP,StateWeights,crit,niters,check) function [StateMu,StateSqrtP,StateWeights] = fit_gaussian_mixture(X,X_weights,StateMu,StateSqrtP,StateWeights,crit,niters,check) % Copyright (C) 2013 Dynare Team % Copyright (C) 2013 Dynare Team % % ... @@ -26,7 +26,7 @@ end ... @@ -26,7 +26,7 @@ end eold = -Inf; eold = -Inf; for n=1:niters for n=1:niters % Calculate posteriors based on old parameters % Calculate posteriors based on old parameters [prior,likelihood,marginal,posterior] = probability(StateMu,StateSqrtP,StateWeights,X); [prior,likelihood,marginal,posterior] = probability3(StateMu,StateSqrtP,StateWeights,X,X_weights); e = sum(log(marginal)); e = sum(log(marginal)); if (n > 1 && abs((e - eold)/eold) < crit) if (n > 1 && abs((e - eold)/eold) < crit) return; return; ... ...
 ... @@ -294,7 +294,7 @@ for t=1:sample_size ... @@ -294,7 +294,7 @@ for t=1:sample_size StateParticles = bsxfun(@plus,StateMuPost(:,i),StateSqrtPPost(:,:,i)*nodes') ; StateParticles = bsxfun(@plus,StateMuPost(:,i),StateSqrtPPost(:,:,i)*nodes') ; IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... StateMuPost,StateSqrtPPost,StateWeightsPost,... StateMuPost,StateSqrtPPost,StateWeightsPost,... StateParticles,H,const_lik,weights,weights_c,ReducedForm,ThreadsOptions) ; StateParticles,H,const_lik,ReducedForm,ThreadsOptions) ; SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights) ; SampleWeights(i) = sum(StateWeightsPost(i)*weights.*IncrementalWeights) ; end end SumSampleWeights = sum(SampleWeights) ; SumSampleWeights = sum(SampleWeights) ; ... @@ -312,13 +312,16 @@ for t=1:sample_size ... @@ -312,13 +312,16 @@ for t=1:sample_size StateParticles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost',number_of_particles) ; StateParticles = importance_sampling(StateMuPost,StateSqrtPPost,StateWeightsPost',number_of_particles) ; IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... IncrementalWeights = gaussian_mixture_densities(Y(:,t),StateMuPrior,StateSqrtPPrior,StateWeightsPrior,... StateMuPost,StateSqrtPPost,StateWeightsPost,... StateMuPost,StateSqrtPPost,StateWeightsPost,... StateParticles,H,const_lik,1/number_of_particles,... StateParticles,H,const_lik,ReducedForm,ThreadsOptions) ; 1/number_of_particles,ReducedForm,ThreadsOptions) ; SampleWeights = IncrementalWeights/number_of_particles ; SampleWeights = IncrementalWeights/number_of_particles ; SumSampleWeights = sum(SampleWeights,1) ; SumSampleWeights = sum(SampleWeights,1) ; SampleWeights = SampleWeights./SumSampleWeights ; SampleWeights = SampleWeights./SumSampleWeights ; lik(t) = log(SumSampleWeights) ; lik(t) = log(SumSampleWeights) ; [StateMu,StateSqrtP,StateWeights] = fit_gaussian_mixture(StateParticles,StateMu,StateSqrtP,StateWeights,0.001,10,1) ; if (ParticleOptions.resampling.status.generic && neff(SampleWeights)
 function [prior,likelihood,C,posterior] = probability3(mu,sqrtP,prior,X,X_weights) % Copyright (C) 2013 Dynare Team % % This file is part of Dynare. % % Dynare is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % (at your option) any later version. % % Dynare is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . [dim,nov] = size(X); M = size(mu,2) ; if nargout>1 likelihood = zeros(M,nov); normfact = (2*pi)^(dim/2); for k=1:M XX = bsxfun(@minus,X,mu(:,k)); S = sqrtP(:,:,k); foo = S \ XX; likelihood(k,:) = exp(-0.5*sum(foo.*foo, 1))/abs((normfact*prod(diag(S)))); end end wlikelihood = bsxfun(@times,X_weights,likelihood) + 1e-99; if nargout>2 C = prior*wlikelihood + 1e-99; end if nargout>3 posterior = bsxfun(@rdivide,bsxfun(@times,prior',wlikelihood),C) + 1e-99 ; posterior = bsxfun(@rdivide,posterior,sum(posterior,1)); end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!