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

Added analytic derivatives for prior distributions.

time-shift
Marco Ratto 2011-11-05 10:32:37 +01:00
parent 8612c7258b
commit 6b2cee2017
5 changed files with 84 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
matlab/lpdfgam.m
View File

@ -1,4 +1,4 @@
function ldens = lpdfgam(x,a,b);
function [ldens,Dldens,D2ldens] = lpdfgam(x,a,b);
% Evaluates the logged GAMMA PDF at x.
%
% INPUTS
@ -37,4 +37,22 @@ if length(a)==1
ldens(idx) = -gammaln(a) - a*log(b) + (a-1)*log(x(idx)) - x(idx)/b ;
else
ldens(idx) = -gammaln(a(idx)) - a(idx).*log(b(idx)) + (a(idx)-1).*log(x(idx)) - x(idx)./b(idx) ;
end
if nargout >1
if length(a)==1
Dldens(idx) = (a-1)./(x(idx)) - ones(length(idx),1)/b ;
else
Dldens(idx) = (a(idx)-1)./(x(idx)) - ones(length(idx),1)./b(idx) ;
end
end
if nargout == 3
if length(a)==1
D2ldens(idx) = -(a-1)./(x(idx)).^2;
else
D2ldens(idx) = -(a(idx)-1)./(x(idx)).^2;
end
end

20
matlab/lpdfgbeta.m
View File

@ -1,4 +1,4 @@
function ldens = lpdfgbeta(x,a,b,aa,bb);
function [ldens,Dldens,D2ldens] = lpdfgbeta(x,a,b,aa,bb);
% Evaluates the logged BETA PDF at x.
%
% INPUTS
@ -38,4 +38,22 @@ if length(a)==1
ldens(idx) = -betaln(a,b) + (a-1)*log(x(idx)-aa) + (b-1)*log(bb-x(idx)) - (a+b-1)*log(bb-aa) ;
else
ldens(idx) = -betaln(a(idx),b(idx)) + (a(idx)-1).*log(x(idx)-aa(idx)) + (b(idx)-1).*log(bb(idx)-x(idx)) - (a(idx)+b(idx)-1).*log(bb(idx)-aa(idx));
end
if nargout >1
if length(a)==1
Dldens(idx) = (a-1)./(x(idx)-aa) - (b-1)./(bb-x(idx)) ;
else
Dldens(idx) = (a(idx)-1)./(x(idx)-aa(idx)) - (b(idx)-1)./(bb(idx)-x(idx));
end
end
if nargout == 3
if length(a)==1
D2ldens(idx) = -(a-1)./(x(idx)-aa).^2 + (b-1)./(bb-x(idx)).^2 ;
else
D2ldens(idx) = -(a(idx)-1)./(x(idx)-aa(idx)).^2 + (b(idx)-1)./(bb(idx)-x(idx)).^2;
end
end

21
matlab/lpdfig1.m
View File

@ -1,4 +1,4 @@
function ldens = lpdfig1(x,s,nu)
function [ldens,Dldens,D2ldens] = lpdfig1(x,s,nu)
% Evaluates the logged INVERSE-GAMMA-1 PDF at x.
%
% X ~ IG1(s,nu) if X = sqrt(Y) where Y ~ IG2(s,nu) and Y = inv(Z) with Z ~ G(nu/2,2/s) (Gamma distribution)
@ -41,4 +41,21 @@ if length(s)==1
ldens(idx) = log(2) - gammaln(.5*nu) - .5*nu*(log(2)-log(s)) - (nu+1)*log(x(idx)) - .5*s./(x(idx).*x(idx)) ;
else
ldens(idx) = log(2) - gammaln(.5*nu(idx)) - .5*nu(idx).*(log(2)-log(s(idx))) - (nu(idx)+1).*log(x(idx)) - .5*s(idx)./(x(idx).*x(idx)) ;
end
end
if nargout >1
if length(s)==1
Dldens(idx) = - (nu+1)./(x(idx)) + s./(x(idx).^3) ;
else
Dldens(idx) = - (nu(idx)+1)./(x(idx)) + s(idx)./(x(idx).^3) ;
end
end
if nargout == 3
if length(s)==1
D2ldens(idx) = (nu+1)./(x(idx).^2) - 3*s.*x(idx).^2./(x(idx).^4) ;
else
D2ldens(idx) = (nu(idx)+1)./(x(idx).^2) - 3*s.*x(idx).^2./(x(idx).^4) ;
end
end

18
matlab/lpdfig2.m
View File

@ -1,4 +1,4 @@
function ldens = lpdfig2(x,s,nu)
function [ldens,Dldens,D2ldens] = lpdfig2(x,s,nu)
% Evaluates the logged INVERSE-GAMMA-2 PDF at x.
%
% X ~ IG2(s,nu) if X = inv(Z) where Z ~ G(nu/2,2/s) (Gamma distribution)
@ -41,4 +41,20 @@ if length(s)==1
ldens(idx) = -gammaln(.5*nu) - (.5*nu)*(log(2)-log(s)) - .5*(nu+2)*log(x(idx)) -.5*s./x(idx);
else
ldens(idx) = -gammaln(.5*nu(idx)) - (.5*nu(idx)).*(log(2)-log(s(idx))) - .5*(nu(idx)+2).*log(x(idx)) -.5*s(idx)./x(idx);
end
if nargout >1
if length(s)==1
Dldens(idx) = - .5*(nu+2)./(x(idx)) + .5*s./x(idx).^2;
else
Dldens(idx) = - .5*(nu(idx)+2)./(x(idx)) + .5*s(idx)./x(idx).^2;
end
end
if nargout == 3
if length(s)==1
D2ldens(idx) = .5*(nu+2)./(x(idx)).^2 - s./x(idx).^3;
else
D2ldens(idx) = .5*(nu(idx)+2)./(x(idx)).^2 - s(idx)./x(idx).^3;
end
end

12
matlab/lpdfnorm.m
View File

@ -1,4 +1,4 @@
function ldens = lpdfnorm(x,a,b)
function [ldens,Dldens,D2ldens] = lpdfnorm(x,a,b)
% Evaluates the logged UNIVARIATE GAUSSIAN PDF at x.
%
% INPUTS
@ -32,4 +32,12 @@ function ldens = lpdfnorm(x,a,b)
if nargin<3, b=1; end
if nargin<2, a=0; end
ldens = -log(b) -.5*log(2*pi) - .5*((x-a)./b).*((x-a)./b) ;
ldens = -log(b) -.5*log(2*pi) - .5*((x-a)./b).*((x-a)./b) ;
if nargout >1
Dldens = - (1/b)*((x-a)/b) ;
end
if nargout == 3
D2ldens = - (1/b)^2 ;
end