Compare commits
2 Commits
master
...
covariance
Author | SHA1 | Date |
---|---|---|
Stéphane Adjemian (Argos) | f3c6328af6 | |
Stéphane Adjemian (Argos) | f1d444e198 |
|
@ -0,0 +1,191 @@
|
||||||
|
function [f, df, d2f, R2] = likelihood_quadratic_approximation(particles, likelihoodvalues)
|
||||||
|
|
||||||
|
% Approximate the shape of the likelihood function with a multivariate second order polynomial.
|
||||||
|
%
|
||||||
|
%
|
||||||
|
% INPUTS
|
||||||
|
% - particles [double] n×p matrix of (p) particles around the estimated posterior mode.
|
||||||
|
% - likelihoodvalues [double] p×1 vector of corresponding values for the likelihood (or posterior kernel).
|
||||||
|
%
|
||||||
|
% OUTPUTS
|
||||||
|
% - f [handle] function handle for the approximated likelihood.
|
||||||
|
% - df [handle] function handle for the gradient of the approximated likelihood.
|
||||||
|
% - d2f [handle] Hessian matrix of the approximated likelihood (constant since we consider a second order multivariate polynomial)
|
||||||
|
% - R2 [double] scalar, goodness of fit measure.
|
||||||
|
%
|
||||||
|
% REMARKS
|
||||||
|
% [1] Function f takes a n×m matrix as input argument (the approximated likelihood is evaluated in m points) and returns a m×1 vector.
|
||||||
|
% [2] Funtion df takes a n×1 vector as input argument (the point where the gradient is computed) and returns a n×1 vector.
|
||||||
|
|
||||||
|
% Copyright © 2024 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 <https://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
|
n = rows(particles); % Number of parmaeters
|
||||||
|
p = columns(particles); % Number of particles
|
||||||
|
q = 1 + n + n*(n+1)/2; % Number of regressors (with a constant)
|
||||||
|
|
||||||
|
if p<=q
|
||||||
|
error('Quadratic approximation requires more than %u particles.', q)
|
||||||
|
end
|
||||||
|
|
||||||
|
%
|
||||||
|
% Build the set of regressors.
|
||||||
|
%
|
||||||
|
|
||||||
|
X = NaN(p, q);
|
||||||
|
X(:,1) = 1; % zero order term
|
||||||
|
X(:,2:n+1) = transpose(particles); % first order terms
|
||||||
|
X(:,n+2:end) = crossproducts(particles); % second order terms
|
||||||
|
|
||||||
|
%
|
||||||
|
% Perform the regression
|
||||||
|
%
|
||||||
|
|
||||||
|
parameters = X\likelihoodvalues(:);
|
||||||
|
|
||||||
|
%
|
||||||
|
% Return a function to evaluate the approximation at x (a n×1 vector).
|
||||||
|
%
|
||||||
|
|
||||||
|
f = @(X) parameters(1) + transpose(X)*parameters(2:n+1) + crossproducts(X)*parameters(n+2:end);
|
||||||
|
|
||||||
|
if nargout>1
|
||||||
|
%
|
||||||
|
% Return a function to evaluate the gradient of the approximation at x (a n×1 vector)
|
||||||
|
%
|
||||||
|
|
||||||
|
df = @(X) parameters(2:n+1) + dcrossproducts(X)*parameters(n+2:end);
|
||||||
|
|
||||||
|
if nargout>2
|
||||||
|
%
|
||||||
|
% Return the hessian matrix of the approximation.
|
||||||
|
%
|
||||||
|
|
||||||
|
d2f = NaN(n,n);
|
||||||
|
|
||||||
|
h = 1;
|
||||||
|
for i=1:n
|
||||||
|
for j=i:n
|
||||||
|
d2f(i,j) = parameters(n+1+h);
|
||||||
|
if ~isequal(j, i)
|
||||||
|
d2f(j,i) = d2f(i,j);
|
||||||
|
end
|
||||||
|
h = h+1;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
if nargout>3
|
||||||
|
%
|
||||||
|
% Return a measure of fit goodness
|
||||||
|
%
|
||||||
|
|
||||||
|
R2 = 1-sum((likelihoodvalues(:)-X*parameters).^2)/sum(demean(likelihoodvalues(:)).^2);
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
function XX = crossproducts(X)
|
||||||
|
% n n
|
||||||
|
% XX*ones(1,(n+1)*n/2) = ∑ xᵢ² + 2 ∑ xᵢxⱼ
|
||||||
|
% i=1 i=1
|
||||||
|
% j>i
|
||||||
|
XX = NaN(columns(X), n*(n+1)/2);
|
||||||
|
column = 1;
|
||||||
|
for i=1:n
|
||||||
|
XX(:,column) = transpose(X(i,:).*X(i,:));
|
||||||
|
column = column+1;
|
||||||
|
for j=i+1:n
|
||||||
|
XX(:,column) = 2*transpose(X(i,:).*X(j,:));
|
||||||
|
column = column+1;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
function xx = dcrossproducts(x)
|
||||||
|
xx = zeros(n, n*(n+1)/2);
|
||||||
|
for i = 1:n
|
||||||
|
base = (i-1)*n-sum(0:i-2);
|
||||||
|
incol = 1;
|
||||||
|
xx(i,base+incol) = 2*x(i);
|
||||||
|
for j = i+1:n
|
||||||
|
incol = incol+1;
|
||||||
|
xx(i,incol) = 2*x(j);
|
||||||
|
end
|
||||||
|
for j=1:i-1
|
||||||
|
base = (j-1)*n-sum(0:j-2)+1;
|
||||||
|
colid = base+i-j;
|
||||||
|
xx(i,colid) = 2*x(j);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
return % --*-- Unit tests --*--
|
||||||
|
|
||||||
|
%@test:1
|
||||||
|
% Create data
|
||||||
|
X = randn(10,1000);
|
||||||
|
Y = 1 + rand(1,10)*X.^2+0.01*randn(1,1000);
|
||||||
|
|
||||||
|
% Perform approximation
|
||||||
|
try
|
||||||
|
[f,df, d2f,R2] = likelihood_quadratic_approximation(X,Y);
|
||||||
|
t(1) = true;
|
||||||
|
catch
|
||||||
|
t(1) = false;
|
||||||
|
end
|
||||||
|
|
||||||
|
% Test returned arguments
|
||||||
|
if t(1)
|
||||||
|
try
|
||||||
|
y = f(randn(10,100));
|
||||||
|
t(2) = true;
|
||||||
|
if ~(rows(y)==100 && columns(y)==1)
|
||||||
|
t(2) = false;
|
||||||
|
end
|
||||||
|
catch
|
||||||
|
t(2) = false;
|
||||||
|
end
|
||||||
|
try
|
||||||
|
dy = df(zeros(10,1));
|
||||||
|
t(3) = true;
|
||||||
|
if ~(rows(dy)==10 && columns(dy)==1)
|
||||||
|
t(3) = false;
|
||||||
|
end
|
||||||
|
catch
|
||||||
|
t(3) = false;
|
||||||
|
end
|
||||||
|
t(4) = true;
|
||||||
|
if ~(rows(d2f)==10 && columns(d2f)==10)
|
||||||
|
t(4) = false
|
||||||
|
end
|
||||||
|
if ~(rows(d2f)==10 && columns(d2f)==10)
|
||||||
|
t(4) = false
|
||||||
|
end
|
||||||
|
t(4) = issymmetric(d2f);
|
||||||
|
t(4) = ispd(d2f);
|
||||||
|
t(5) = isscalar(R2);
|
||||||
|
t(5) = (R2>0) & (R2<1); % Note that in a nonlinear model nothing ensures that these inequalities are satisfied.
|
||||||
|
end
|
||||||
|
|
||||||
|
T = all(t);
|
||||||
|
%@eof:1
|
||||||
|
|
||||||
|
end
|
|
@ -1,30 +1,15 @@
|
||||||
function [test, penalty] = ispd(A)
|
function [test, penalty] = ispd(A)
|
||||||
|
|
||||||
%@info:
|
% Test if the square matrix A is positive definite.
|
||||||
%! @deftypefn {Function File} {[@var{test}, @var{penalty} =} ispd (@var{A})
|
%
|
||||||
%! @anchor{ispd}
|
% INPUTS
|
||||||
%! @sp 1
|
% - A [double] n×n matrix.
|
||||||
%! Tests if the square matrix @var{A} is positive definite.
|
%
|
||||||
%! @sp 2
|
% OUTPUTS
|
||||||
%! @strong{Inputs}
|
% - test [logical] scalar, true if and only if matrix A is positive definite (and symmetric...)
|
||||||
%! @sp 1
|
% - penalty [double] scalar, absolute value of the uum of the negative eigenvalues of A. This output argument is optional.
|
||||||
%! @table @ @var
|
|
||||||
%! @item A
|
|
||||||
%! A square matrix.
|
|
||||||
%! @end table
|
|
||||||
%! @sp 2
|
|
||||||
%! @strong{Outputs}
|
|
||||||
%! @sp 1
|
|
||||||
%! @table @ @var
|
|
||||||
%! @item test
|
|
||||||
%! Integer scalar equal to 1 if @var{A} is a positive definite sqquare matrix, 0 otherwise.
|
|
||||||
%! @item penalty
|
|
||||||
%! Absolute value of the uum of the negative eigenvalues of A. This output argument is optional.
|
|
||||||
%! @end table
|
|
||||||
%! @end deftypefn
|
|
||||||
%@eod:
|
|
||||||
|
|
||||||
% Copyright © 2007-2022 Dynare Team
|
% Copyright © 2007-2024 Dynare Team
|
||||||
%
|
%
|
||||||
% This file is part of Dynare.
|
% This file is part of Dynare.
|
||||||
%
|
%
|
||||||
|
@ -54,6 +39,7 @@ if nargout>1
|
||||||
if isoctave && any(any(~isfinite(A))) % workaround for https://savannah.gnu.org/bugs/index.php?63082
|
if isoctave && any(any(~isfinite(A))) % workaround for https://savannah.gnu.org/bugs/index.php?63082
|
||||||
penalty = 1;
|
penalty = 1;
|
||||||
else
|
else
|
||||||
|
% TODO: the penalty is only concerned with negative eigenvalues, we should also consider the case of non symmetric matrix A.
|
||||||
a = diag(eig(A));
|
a = diag(eig(A));
|
||||||
k = find(a<0);
|
k = find(a<0);
|
||||||
if k > 0
|
if k > 0
|
||||||
|
@ -61,4 +47,4 @@ if nargout>1
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
Loading…
Reference in New Issue