probability.m 1.28 KB
 Stéphane Adjemian committed Nov 10, 2014 1 2 ``````function [prior,likelihood,C,posterior] = probability(mu,sqrtP,prior,X) `````` Stéphane Adjemian committed May 18, 2017 3 ``````% Copyright (C) 2013-2017 Dynare Team `````` Stéphane Adjemian committed Nov 10, 2014 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ``````% % 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 . `````` Stéphane Adjemian committed May 18, 2017 20 ``````[dim,nov] = size(X); `````` Stéphane Adjemian committed Nov 10, 2014 21 22 ``````M = size(mu,2) ; if nargout>1 `````` Stéphane Adjemian committed May 18, 2017 23 24 25 26 27 28 29 30 `````` 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 `````` Stéphane Adjemian committed Nov 10, 2014 31 32 33 ``````end likelihood = likelihood + 1e-99; if nargout>2 `````` Stéphane Adjemian committed May 18, 2017 34 `````` C = prior*likelihood + 1e-99; `````` Stéphane Adjemian committed Nov 10, 2014 35 36 ``````end if nargout>3 `````` Stéphane Adjemian committed May 18, 2017 37 38 `````` posterior = bsxfun(@rdivide,bsxfun(@times,prior',likelihood),C) + 1e-99 ; posterior = bsxfun(@rdivide,posterior,sum(posterior,1)); `````` Stéphane Adjemian committed Nov 10, 2014 39 ``end``