v3+4: mr_hessian.m corrected by M. Ratto
git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@592 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
parent
c642044b39
commit
e62b80c658
|
@ -1,33 +1,68 @@
|
||||||
% Copyright (C) 2001 Michel Juillard
|
% Copyright (C) 2004 Marco Ratto
|
||||||
|
% adapted from Michel Juillard original rutine hessian.m
|
||||||
%
|
%
|
||||||
% computes second order partial derivatives
|
% [hessian_mat, gg, htol1, ihh, hh_mat0] = mr_hessian(func,x,hflag,htol0,varargin)
|
||||||
% uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27 p. 884
|
|
||||||
%
|
%
|
||||||
% Adapted by M. Ratto from original M. Juillard routine
|
% numerical gradient and Hessian, with 'automatic' check of numerical
|
||||||
|
% error
|
||||||
|
%
|
||||||
|
% func = name of the function: func must give two outputs:
|
||||||
|
% - the log-likelihood AND the single contributions at times t=1,...,T
|
||||||
|
% of the log-likelihood to compute outer product gradient
|
||||||
|
% x = parameter values
|
||||||
|
% hflag = 0, Hessian computed with outer product gradient, one point
|
||||||
|
% increments for partial derivatives in gradients
|
||||||
|
% hflag = 1, 'mixed' Hessian: diagonal elements computed with numerical second order derivatives
|
||||||
|
% with correlation structure as from outer product gradient;
|
||||||
|
% two point evaluation of derivatives for partial derivatives
|
||||||
|
% in gradients
|
||||||
|
% hflag = 2, full numerical Hessian, computes second order partial derivatives
|
||||||
|
% uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27
|
||||||
|
% p. 884.
|
||||||
|
% htol0 = 'precision' of increment of function values for numerical
|
||||||
|
% derivatives
|
||||||
|
%
|
||||||
|
% varargin: other parameters of func
|
||||||
%
|
%
|
||||||
function [hessian_mat, gg] = mr_hessian(func,x,hflag,varargin)
|
|
||||||
global gstep_
|
|
||||||
persistent h1
|
|
||||||
|
|
||||||
|
function [hessian_mat, gg, htol1, ihh, hh_mat0] = mr_hessian(func,x,hflag,htol0,varargin)
|
||||||
|
global gstep_ bayestopt_
|
||||||
|
persistent h1 htol
|
||||||
|
|
||||||
|
if isempty(htol), htol = 1.e-4; end
|
||||||
func = str2func(func);
|
func = str2func(func);
|
||||||
f0=feval(func,x,varargin{:});
|
[f0, ff0]=feval(func,x,varargin{:});
|
||||||
n=size(x,1);
|
n=size(x,1);
|
||||||
|
h2=bayestopt_.ub-bayestopt_.lb;
|
||||||
%h1=max(abs(x),gstep_*ones(n,1))*eps^(1/3);
|
%h1=max(abs(x),gstep_*ones(n,1))*eps^(1/3);
|
||||||
%h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/6);
|
%h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/6);
|
||||||
if isempty(h1),
|
if isempty(h1),
|
||||||
h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/4);
|
h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/4);
|
||||||
end
|
end
|
||||||
htol = 1.e-4;
|
|
||||||
|
if htol0<htol,
|
||||||
|
htol=htol0;
|
||||||
|
end
|
||||||
xh1=x;
|
xh1=x;
|
||||||
f1=zeros(size(f0,1),n);
|
f1=zeros(size(f0,1),n);
|
||||||
f_1=f1;
|
f_1=f1;
|
||||||
|
ff1=zeros(size(ff0));
|
||||||
|
ff_1=ff1;
|
||||||
|
|
||||||
for i=1:n,
|
%for i=1:n,
|
||||||
|
i=0;
|
||||||
|
while i<n,
|
||||||
|
i=i+1;
|
||||||
|
h10=h1(i);
|
||||||
|
hcheck=0;
|
||||||
|
dx=[];
|
||||||
xh1(i)=x(i)+h1(i);
|
xh1(i)=x(i)+h1(i);
|
||||||
fx=feval(func,xh1,varargin{:});
|
[fx, ffx]=feval(func,xh1,varargin{:});
|
||||||
if abs(fx-f0)<htol | abs(fx-f0)>(5*htol),
|
it=1;
|
||||||
|
dx=(fx-f0);
|
||||||
|
ic=0;
|
||||||
|
if abs(dx)>(2*htol),
|
||||||
c=mr_nlincon(xh1,varargin{:});
|
c=mr_nlincon(xh1,varargin{:});
|
||||||
ic=0;
|
|
||||||
while c
|
while c
|
||||||
h1(i)=h1(i)*0.9;
|
h1(i)=h1(i)*0.9;
|
||||||
xh1(i)=x(i)+h1(i);
|
xh1(i)=x(i)+h1(i);
|
||||||
|
@ -35,53 +70,111 @@ for i=1:n,
|
||||||
ic=1;
|
ic=1;
|
||||||
end
|
end
|
||||||
if ic,
|
if ic,
|
||||||
fx=feval(func,xh1,varargin{:});
|
[fx, ffx]=feval(func,xh1,varargin{:});
|
||||||
|
dx=(fx-f0);
|
||||||
end
|
end
|
||||||
|
end
|
||||||
icount = 0;
|
|
||||||
while abs(fx-f0)<htol & icount< 10,
|
icount = 0;
|
||||||
icount=icount+1;
|
h0=h1(i);
|
||||||
h1(i)=min(0.3*abs(x(i)), 1.e-3/(abs(fx-f0)/h1(i)));
|
while (abs(dx(it))<0.5*htol | abs(dx(it))>(2*htol)) & icount<10 & ic==0,
|
||||||
|
%while abs(dx(it))<0.5*htol & icount< 10 & ic==0,
|
||||||
|
icount=icount+1;
|
||||||
|
%if abs(dx(it)) ~= 0,
|
||||||
|
if abs(dx(it))<0.5*htol
|
||||||
|
if abs(dx(it)) ~= 0,
|
||||||
|
h1(i)=min(0.3*abs(x(i)), 0.9*htol/abs(dx(it))*h1(i));
|
||||||
|
else
|
||||||
|
h1(i)=2.1*h1(i);
|
||||||
|
end
|
||||||
xh1(i)=x(i)+h1(i);
|
xh1(i)=x(i)+h1(i);
|
||||||
c=mr_nlincon(xh1,varargin{:});
|
c=mr_nlincon(xh1,varargin{:});
|
||||||
while c
|
while c
|
||||||
h1(i)=h1(i)*0.9;
|
h1(i)=h1(i)*0.9;
|
||||||
xh1(i)=x(i)+h1(i);
|
xh1(i)=x(i)+h1(i);
|
||||||
c=mr_nlincon(xh1,varargin{:});
|
c=mr_nlincon(xh1,varargin{:});
|
||||||
end
|
ic=1;
|
||||||
fx=feval(func,xh1,varargin{:});
|
end
|
||||||
|
[fx, ffx]=feval(func,xh1,varargin{:});
|
||||||
end
|
end
|
||||||
while abs(fx-f0)>(5*htol),
|
if abs(dx(it))>(2*htol),
|
||||||
h1(i)=h1(i)*0.5;
|
h1(i)= htol/abs(dx(it))*h1(i);
|
||||||
xh1(i)=x(i)+h1(i);
|
xh1(i)=x(i)+h1(i);
|
||||||
fx=feval(func,xh1,varargin{:});
|
[fx, ffx]=feval(func,xh1,varargin{:});
|
||||||
|
while (fx-f0)==0,
|
||||||
|
h1(i)= h1(i)*2;
|
||||||
|
xh1(i)=x(i)+h1(i);
|
||||||
|
[fx, ffx]=feval(func,xh1,varargin{:});
|
||||||
|
ic=1;
|
||||||
|
end
|
||||||
end
|
end
|
||||||
|
it=it+1;
|
||||||
|
dx(it)=(fx-f0);
|
||||||
|
h0(it)=h1(i);
|
||||||
|
if h1(i)<1.e-12*min(1,h2(i)),
|
||||||
|
ic=1;
|
||||||
|
hcheck=1;
|
||||||
|
end
|
||||||
|
%else
|
||||||
|
% h1(i)=1;
|
||||||
|
% ic=1;
|
||||||
|
%end
|
||||||
end
|
end
|
||||||
|
% if (it>2 & dx(1)<10^(log10(htol)/2)) ,
|
||||||
|
% [dum, is]=sort(h0);
|
||||||
|
% if find(diff(sign(diff(dx(is)))));
|
||||||
|
% hcheck=1;
|
||||||
|
% end
|
||||||
|
% elseif (it>3 & dx(1)>10^(log10(htol)/2)) ,
|
||||||
|
% [dum, is]=sort(h0);
|
||||||
|
% if find(diff(sign(diff(dx(is(1:end-1))))));
|
||||||
|
% hcheck=1;
|
||||||
|
% end
|
||||||
|
% end
|
||||||
f1(:,i)=fx;
|
f1(:,i)=fx;
|
||||||
xh1(i)=x(i)-h1(i);
|
ff1=ffx;
|
||||||
c=mr_nlincon(xh1,varargin{:});
|
if hflag, % two point based derivatives
|
||||||
ic=0;
|
|
||||||
while c
|
|
||||||
h1(i)=h1(i)*0.9;
|
|
||||||
xh1(i)=x(i)-h1(i);
|
xh1(i)=x(i)-h1(i);
|
||||||
c=mr_nlincon(xh1,varargin{:});
|
c=mr_nlincon(xh1,varargin{:});
|
||||||
ic = 1;
|
ic=0;
|
||||||
end
|
while c
|
||||||
fx=feval(func,xh1,varargin{:});
|
h1(i)=h1(i)*0.9;
|
||||||
f_1(:,i)=fx;
|
xh1(i)=x(i)-h1(i);
|
||||||
if ic,
|
c=mr_nlincon(xh1,varargin{:});
|
||||||
xh1(i)=x(i)+h1(i);
|
ic = 1;
|
||||||
f1(:,i)=feval(func,xh1,varargin{:});
|
end
|
||||||
|
[fx, ffx]=feval(func,xh1,varargin{:});
|
||||||
|
f_1(:,i)=fx;
|
||||||
|
ff_1=ffx;
|
||||||
|
if ic,
|
||||||
|
xh1(i)=x(i)+h1(i);
|
||||||
|
[f1(:,i), ff1]=feval(func,xh1,varargin{:});
|
||||||
|
end
|
||||||
|
ggh(:,i)=(ff1-ff_1)./(2.*h1(i));
|
||||||
|
else
|
||||||
|
ggh(:,i)=(ff1-ff0)./h1(i);
|
||||||
end
|
end
|
||||||
xh1(i)=x(i);
|
xh1(i)=x(i);
|
||||||
|
if hcheck & htol<1,
|
||||||
|
htol=min(1,max(min(abs(dx))*2,htol*10));
|
||||||
|
h1(i)=h10;
|
||||||
|
i=0;
|
||||||
|
end
|
||||||
|
save hess
|
||||||
end
|
end
|
||||||
|
|
||||||
h_1=h1;
|
h_1=h1;
|
||||||
xh1=x;
|
xh1=x;
|
||||||
xh_1=xh1;
|
xh_1=xh1;
|
||||||
gg=(f1'-f_1')./(2.*h1);
|
|
||||||
|
|
||||||
if hflag,
|
if hflag,
|
||||||
|
gg=(f1'-f_1')./(2.*h1);
|
||||||
|
else
|
||||||
|
gg=(f1'-f0)./h1;
|
||||||
|
end
|
||||||
|
|
||||||
|
if hflag==2,
|
||||||
|
gg=(f1'-f_1')./(2.*h1);
|
||||||
hessian_mat = zeros(size(f0,1),n*n);
|
hessian_mat = zeros(size(f0,1),n*n);
|
||||||
for i=1:n
|
for i=1:n
|
||||||
if i > 1
|
if i > 1
|
||||||
|
@ -96,35 +189,80 @@ if hflag,
|
||||||
xh_1(i)=x(i)-h1(i);
|
xh_1(i)=x(i)-h1(i);
|
||||||
xh_1(j)=x(j)-h_1(j);
|
xh_1(j)=x(j)-h_1(j);
|
||||||
%hessian_mat(:,(i-1)*n+j)=-(-feval(func,xh1,varargin{:})-feval(func,xh_1,varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
|
%hessian_mat(:,(i-1)*n+j)=-(-feval(func,xh1,varargin{:})-feval(func,xh_1,varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
|
||||||
temp1 = feval(func,xh1,varargin{:});
|
%temp1 = feval(func,xh1,varargin{:});
|
||||||
%c=mr_nlincon(xh1,varargin{:});
|
c=mr_nlincon(xh1,varargin{:});
|
||||||
%if c, disp( ['hessian warning cross ', num2str(c) ]), end
|
lam=1;
|
||||||
|
while c,
|
||||||
|
lam=lam*0.9;
|
||||||
|
xh1(i)=x(i)+h1(i)*lam;
|
||||||
|
xh1(j)=x(j)+h_1(j)*lam;
|
||||||
|
%disp( ['hessian warning cross ', num2str(c) ]),
|
||||||
|
c=mr_nlincon(xh1,varargin{:});
|
||||||
|
end
|
||||||
|
temp1 = f0+(feval(func,xh1,varargin{:})-f0)/lam;
|
||||||
|
|
||||||
|
%temp2 = feval(func,xh_1,varargin{:});
|
||||||
|
c=mr_nlincon(xh_1,varargin{:});
|
||||||
|
while c,
|
||||||
|
lam=lam*0.9;
|
||||||
|
xh_1(i)=x(i)-h1(i)*lam;
|
||||||
|
xh_1(j)=x(j)-h_1(j)*lam;
|
||||||
|
%disp( ['hessian warning cross ', num2str(c) ]),
|
||||||
|
c=mr_nlincon(xh_1,varargin{:});
|
||||||
|
end
|
||||||
|
temp2 = f0+(feval(func,xh_1,varargin{:})-f0)/lam;
|
||||||
|
|
||||||
temp2 = feval(func,xh_1,varargin{:});
|
|
||||||
%c=mr_nlincon(xh_1,varargin{:});
|
|
||||||
%if c, disp( ['hessian warning cross ', num2str(c) ]), end
|
|
||||||
hessian_mat(:,(i-1)*n+j)=-(-temp1 -temp2+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
|
hessian_mat(:,(i-1)*n+j)=-(-temp1 -temp2+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
|
||||||
xh1(i)=x(i);
|
xh1(i)=x(i);
|
||||||
xh1(j)=x(j);
|
xh1(j)=x(j);
|
||||||
xh_1(i)=x(i);
|
xh_1(i)=x(i);
|
||||||
xh_1(j)=x(j);
|
xh_1(j)=x(j);
|
||||||
j=j+1;
|
j=j+1;
|
||||||
|
save hess
|
||||||
end
|
end
|
||||||
i=i+1;
|
i=i+1;
|
||||||
end
|
end
|
||||||
|
|
||||||
else
|
elseif hflag==1,
|
||||||
hessian_mat = zeros(size(f0,1),n*n);
|
hessian_mat = zeros(size(f0,1),n*n);
|
||||||
for i=1:n,
|
for i=1:n,
|
||||||
dum = (f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
|
dum = (f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
|
||||||
if dum>0,
|
if dum>eps,
|
||||||
hessian_mat(:,(i-1)*n+i)=dum;
|
hessian_mat(:,(i-1)*n+i)=dum;
|
||||||
else
|
else
|
||||||
hessian_mat(:,(i-1)*n+i)=gg(i)^2;
|
hessian_mat(:,(i-1)*n+i)=max(eps, gg(i)^2);
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
%hessian_mat2=hh_mat(:)';
|
||||||
|
end
|
||||||
|
|
||||||
|
gga=ggh.*kron(ones(size(ff1)),2.*h1'); % re-scaled gradient
|
||||||
|
hh_mat=gga'*gga; % rescaled outer product hessian
|
||||||
|
hh_mat0=ggh'*ggh; % outer product hessian
|
||||||
|
A=diag(2.*h1); % rescaling matrix
|
||||||
|
if hflag>0 & min(eig(reshape(hessian_mat,n,n)))>0,
|
||||||
|
hh0 = A*reshape(hessian_mat,n,n)*A'; %rescaled second order derivatives
|
||||||
|
sd0=sqrt(diag(hh0)); %rescaled 'standard errors' using second order derivatives
|
||||||
|
sd=sqrt(diag(hh_mat)); %rescaled 'standard errors' using outer product
|
||||||
|
hh_mat=hh_mat./(sd*sd').*(sd0*sd0'); %rescaled outer product with 'true' std's
|
||||||
|
hh0 = reshape(hessian_mat,n,n); % second order derivatives
|
||||||
|
sd0=sqrt(diag(hh0)); % 'standard errors' using second order derivatives
|
||||||
|
sd=sqrt(diag(hh_mat0)); % 'standard errors' using outer product
|
||||||
|
hh_mat0=hh_mat0./(sd*sd').*(sd0*sd0'); % rescaled outer product with 'true' std's
|
||||||
|
end
|
||||||
|
igg=inv(hh_mat); % inverted rescaled outer product hessian
|
||||||
|
ihh=A'*igg*A; % inverted outer product hessian
|
||||||
|
if hflag==0,
|
||||||
|
hessian_mat=hh_mat0(:);
|
||||||
|
end
|
||||||
|
|
||||||
|
if isnan(hessian_mat),
|
||||||
|
hh_mat0=eye(length(hh_mat0));
|
||||||
|
ihh=hh_mat0;
|
||||||
|
hessian_mat=hh_mat0(:);
|
||||||
end
|
end
|
||||||
hh1=h1;
|
hh1=h1;
|
||||||
|
htol1=htol;
|
||||||
save hess
|
save hess
|
||||||
% 11/25/03 SA Created from Hessian_sparse (removed sparse)
|
% 11/25/03 SA Created from Hessian_sparse (removed sparse)
|
||||||
|
|
||||||
|
|
Loading…
Reference in New Issue