diff --git a/matlab/@dprior/bounds.m b/matlab/@dprior/bounds.m
new file mode 100644
index 0000000000000000000000000000000000000000..d4047ba6e36b64d1641fe0d1189fe9844ab82879
--- /dev/null
+++ b/matlab/@dprior/bounds.m
@@ -0,0 +1,225 @@
+function o = bounds(o, truncation, inplace)
+
+% Return true iff d is an admissible draw in a distribution characterized by o.
+%
+% INPUTS
+% - o [dprior] Distribution specification for a n×1 vector of independent continuous random variables
+% - truncation [double] scalar, prior mass excluded in left and right parts of the prior distribution. Defaut: 0 if inplace==true.
+% - inplace [logical] scalar, if true modifies o in place.
+%
+% OUTPUTS
+% - lub [double] n×2 matrix, lower and upper bounds of the prior support
+
+% 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/>.
+
+if nargin<3 || isempty(inplace)
+ inplace = false;
+end
+
+if inplace
+ if nargin<2 || isempty(truncation)
+ truncation = .0;
+ end
+else
+ if ~isempty(o.lb) && (nargin<2 || isempty(truncation) || abs(truncation-o.trunc)<eps(1))
+ lub = [o.lb, o.ub];
+ o = lub;
+ return
+ end
+end
+
+assert(truncation>=0 && truncation<=1, 'Second input argument must be non negative and not larger than one.')
+
+lub = zeros(length(o.p6), 2);
+
+if o.isbeta
+ if truncation<eps(1)
+ lub(o.idbeta,1) = o.p3(o.idbeta);
+ lub(o.idbeta,2) = o.p4(o.idbeta);
+ else
+ lub(o.idbeta,1) = betainv(truncation, o.p6(o.idbeta), o.p7(o.idbeta)).*(o.p4(o.idbeta)-o.p3(o.idbeta))+o.p3(o.idbeta);
+ lub(o.idbeta,2) = betainv(1.0-truncation, o.p6(o.idbeta), o.p7(o.idbeta)).*(o.p4(o.idbeta)-o.p3(o.idbeta))+o.p3(o.idbeta);
+ end
+end
+
+if o.isgamma
+ if truncation<eps(1)
+ lub(o.idgamma,1) = o.p3(o.idgamma);
+ lub(o.idgamma,2) = Inf;
+ else
+ lub(o.idgamma,1) = gaminv(truncation, o.p6(o.idgamma), o.p7(o.idgamma))+o.p3(o.idgamma);
+ lub(o.idgamma,2) = gaminv(1.0-truncation, o.p6(o.idgamma), o.p7(o.idgamma))+o.p3(o.idgamma);
+ end
+end
+
+if o.isgaussian
+ if truncation<eps(1)
+ lub(o.idgaussian,1) = max(-Inf, o.p3(o.idgaussian));
+ lub(o.idgaussian,2) = min(Inf, o.p4(o.idgaussian));
+ else
+ lub(o.idgaussian,1) = max(norminv(truncation, o.p6(o.idgaussian), o.p7(o.idgaussian)), o.p3(o.idgaussian));
+ lub(o.idgaussian,2) = min(norminv(1-truncation, o.p6(o.idgaussian), o.p7(o.idgaussian)), o.p4(o.idgaussian));
+ end
+end
+
+if o.isinvgamma1
+ if truncation<eps(1)
+ lub(o.idinvgamma1,1) = o.p3(o.idinvgamma1);
+ lub(o.idinvgamma1,2) = Inf;
+ else
+ lub(o.idinvgamma1,1) = 1.0./sqrt(gaminv(1.0-truncation, o.p7(o.idinvgamma1)/2.0, 2.0./o.p6(o.idinvgamma1)))+o.p3(o.idinvgamma1);
+ lub(o.idinvgamma1,2) = 1.0./sqrt(gaminv(truncation, o.p7(o.idinvgamma1)/2.0, 2.0./o.p6(o.idinvgamma1)))+o.p3(o.idinvgamma1);
+ end
+end
+
+if o.isuniform
+ if truncation<eps(1)
+ lub(o.iduniform,1) = o.p6(o.iduniform);
+ lub(o.iduniform,2) = o.p7(o.iduniform);
+ else
+ lub(o.iduniform,1) = o.p6(o.iduniform)+(o.p7(o.iduniform)-o.p6(o.iduniform))*truncation;
+ lub(o.iduniform,2) = o.p7(o.iduniform)-(o.p7(o.iduniform)-o.p6(o.iduniform))*truncation;
+ end
+end
+
+if o.isinvgamma2
+ if truncation<eps(1)
+ lub(o.idinvgamma2,1) = o.p3(o.idinvgamma2);
+ lub(o.idinvgamma2,2) = Inf;
+ else
+ lub(o.idinvgamma2,1) = 1.0./gaminv(1.0-truncation, o.p7(o.idinvgamma2)/2.0, 2.0./o.p6(o.idinvgamma2))+o.p3(o.idinvgamma2);
+ lub(o.idinvgamma2,2) = 1.0./gaminv(truncation, o.p7(o.idinvgamma2)/2.0, 2.0./o.p6(o.idinvgamma2))+o.p3(o.idinvgamma2);
+ end
+end
+
+if o.isweibull
+ if truncation<eps(1)
+ lub(o.idweibull,1) = o.p3(o.idweibull);
+ lub(o.idweibull,2) = Inf;
+ else
+ lub(o.idweibull,1) = o.p3(o.idweibull)+wblinv(truncation, o.p6(o.idweibull), o.p7(o.idweibull));
+ lub(o.idweibull,2) = o.p3(o.idweibull)+wblinv(1.0-truncation, o.p6(o.idweibull), o.p7(o.idweibull));
+ end
+end
+
+if inplace
+ o.lb = lub(:,1);
+ o.ub = lub(:,2);
+else
+ o = lub;
+end
+
+return % --*-- Unit tests --*--
+
+%@test:1
+% Fill global structures with required fields...
+prior_trunc = 1e-10;
+p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1); % Prior shape
+p1 = .4*ones(14,1); % Prior mean
+p2 = .2*ones(14,1); % Prior std.
+p3 = NaN(14,1);
+p4 = NaN(14,1);
+p5 = NaN(14,1);
+p6 = NaN(14,1);
+p7 = NaN(14,1);
+
+for i=1:14
+ switch p0(i)
+ case 1
+ % Beta distribution
+ p3(i) = 0;
+ p4(i) = 1;
+ [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
+ p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
+ case 2
+ % Gamma distribution
+ p3(i) = 0;
+ p4(i) = Inf;
+ [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
+ p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
+ case 3
+ % Normal distribution
+ p3(i) = -Inf;
+ p4(i) = Inf;
+ p6(i) = p1(i);
+ p7(i) = p2(i);
+ p5(i) = p1(i);
+ case 4
+ % Inverse Gamma (type I) distribution
+ p3(i) = 0;
+ p4(i) = Inf;
+ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
+ p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
+ case 5
+ % Uniform distribution
+ [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
+ p3(i) = p6(i);
+ p4(i) = p7(i);
+ p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
+ case 6
+ % Inverse Gamma (type II) distribution
+ p3(i) = 0;
+ p4(i) = Inf;
+ [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
+ p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
+ case 8
+ % Weibull distribution
+ p3(i) = 0;
+ p4(i) = Inf;
+ [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
+ p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
+ otherwise
+ error('This density is not implemented!')
+ end
+end
+
+BayesInfo.pshape = p0;
+BayesInfo.p1 = p1;
+BayesInfo.p2 = p2;
+BayesInfo.p3 = p3;
+BayesInfo.p4 = p4;
+BayesInfo.p5 = p5;
+BayesInfo.p6 = p6;
+BayesInfo.p7 = p7;
+
+% Call the tested routine
+try
+ % Instantiate dprior object
+ o = dprior(BayesInfo, prior_trunc, false);
+ % Save lower and upper bounds
+ lb = o.lb;
+ ub = o.ub;
+ % Call the tested method
+ lub = o.bounds(.1, false);
+ t(1) = true;
+catch
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = all(abs(o.lb-lb)<2*eps(1));
+ t(3) = all(abs(o.ub-ub)<2*eps(1));
+ t(4) = all(lub(:,1)>=lb);
+ t(5) = all(lub(:,2)<=ub);
+end
+
+T = all(t);
+
+%@eof:1
+
+end
diff --git a/matlab/@dprior/dprior.m b/matlab/@dprior/dprior.m
index 0f81c639b2a69468af8e2b6b6c3ac68b95ded6f0..4978447e280e39c312c658e8817618a27c928d4c 100644
--- a/matlab/@dprior/dprior.m
+++ b/matlab/@dprior/dprior.m
@@ -28,7 +28,8 @@ classdef dprior < handle
p11 = []; % Prior median
lb = []; % Truncated prior lower bound.
ub = []; % Truncated prior upper bound.
- name = {}; % Name of the parameter
+ trunc = .0; % Truncated mass.
+ name = {}; % Name of the parameter.
iduniform = []; % Index for the uniform priors.
idgaussian = []; % Index for the gaussian priors.
idgamma = []; % Index for the gamma priors.
@@ -52,11 +53,11 @@ classdef dprior < handle
%
% INPUTS
% - bayestopt_ [struct] Informations about the prior distribution, aka bayestopt_.
- % - PriorTrunc [double] scalar, probability mass to be excluded, aka options_.prior_trunc
- % - Uniform [logical] scalar, produce uniform random deviates on the prior support.
+ % - PriorTrunc [double] scalar, probability mass to be excluded, aka options_.prior_trunc
+ % - Uniform [logical] scalar, produce uniform random deviates on the prior support.
%
% OUTPUTS
- % - o [dprior] scalar, prior object.
+ % - o [dprior] scalar, prior object.
%
% REQUIREMENTS
% None.
@@ -71,9 +72,6 @@ classdef dprior < handle
if isfield(bayestopt_, 'p6'), o.p6 = bayestopt_.p6; end
if isfield(bayestopt_, 'p7'), o.p7 = bayestopt_.p7; end
if isfield(bayestopt_, 'p11'), o.p11 = bayestopt_.p11; end
- bounds = prior_bounds(bayestopt_, PriorTrunc);
- o.lb = bounds.lb;
- o.ub = bounds.ub;
if nargin>2 && Uniform
prior_shape = repmat(5, length(o.p6), 1);
else
@@ -107,6 +105,8 @@ classdef dprior < handle
if ~isempty(o.idweibull)
o.isweibull = true;
end
+ o.trunc = PriorTrunc;
+ o.bounds(o.trunc, true);
end % dprior (constructor)
end % methods
diff --git a/matlab/@dprior/subsref.m b/matlab/@dprior/subsref.m
index faea66d53ba0399125922eab70b73132368dd9c5..692f813d0055dafeeebe58a5e6ac0a67e86f3594 100644
--- a/matlab/@dprior/subsref.m
+++ b/matlab/@dprior/subsref.m
@@ -21,11 +21,13 @@ function p = subsref(o, S)
switch S(1).type
case '.'
- if ismember(S(1).subs, {'p1','p2','p3','p4','p5','p6','p7','lb','ub'})
+ if ismember(S(1).subs, {'p1','p2','p3','p4','p5','p6','p7','lb','ub','trunc'})
p = builtin('subsref', o, S(1));
+ elseif isequal(S(1).subs, 'bounds') && length(S)==1
+ p = bounds(o, [], false);
elseif ismember(S(1).subs, {'draw','length'})
p = feval(S(1).subs, o);
- elseif ismember(S(1).subs, {'draws', 'density', 'densities', 'moments', 'admissible'})
+ elseif ismember(S(1).subs, {'draws', 'density', 'densities', 'moments', 'admissible', 'bounds'})
p = feval(S(1).subs, o , S(2).subs{:});
elseif ismember(S(1).subs, {'mean', 'median', 'variance', 'mode'})
if (length(S)==2 && isempty(S(2).subs)) || length(S)==1
diff --git a/matlab/missing/stats/wblinv.m b/matlab/missing/stats/wblinv.m
index aee1f1706ea5f505a68323bd1c2dcddee0fff30c..bc4024df14dff1efa7ee3fd7f546f930c72de43f 100644
--- a/matlab/missing/stats/wblinv.m
+++ b/matlab/missing/stats/wblinv.m
@@ -27,43 +27,69 @@ function t = wblinv(proba, scale, shape)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
+%
% Check input arguments.
+%
if nargin<3
error('Three input arguments required!')
end
-if ~isnumeric(proba) || ~isscalar(proba) || ~isreal(proba) || proba<0 || proba>1
- error('First input argument must be a real scalar between 0 and 1 (probability)!')
+if ~isnumeric(proba) || ~(isscalar(proba) || isvector(proba)) || ~isreal(proba) || any(proba<0) || any(proba>1)
+ error('First input argument must be a real scalar or vector between 0 and 1 (probability)!')
end
-if ~isnumeric(scale) || ~isscalar(scale) || ~isreal(scale) || scale<=0
- error('Second input argument must be a real positive scalar (scale parameter of the Weibull distribution)!')
+if ~isnumeric(scale) || ~(isscalar(scale) || isvector(scale)) || ~isreal(scale) || any(scale<=0)
+ error('Second input argument must be a real positive scalar or a vector (scale parameter of the Weibull distribution)!')
end
-if ~isnumeric(shape) || ~isscalar(shape) || ~isreal(shape) || shape<=0
- error('Third input argument must be a real positive scalar (shape parameter of the Weibull distribution)!')
+if ~isnumeric(shape) || ~(isscalar(shape) || isvector(shape)) || ~isreal(shape) || any(shape<=0)
+ error('Third input argument must be a real positive scalar or vector (shape parameter of the Weibull distribution)!')
end
+if isvector(shape) && isvector(scale) && ~(isscalar(shape) || isscalar(scale))
+ if length(shape)~=length(scale)
+ error('Second and third arguments (scale and shape) must be vectors of same lengths.')
+ end
+elseif isscalar(scale) && isvector(shape) && ~isscalar(shape)
+ scale = repmat(scale, length(shape), 1);
+elseif isscalar(shape) && isvector(scale) && ~isscalar(scale)
+ shape = repmat(shape, length(scale), 1);
+elseif isscalar(shape) && isscalar(scale)
+ if isvector(proba) && ~isscalar(proba)
+ scale = repmat(scale, length(proba), 1);
+ shape = repmat(shape, length(proba), 1);
+ end
+else
+ error('Second and third arguments (scale and shape) must be vectors or scalars.')
+end
-if proba<2*eps()
- t = 0;
- return
+if isscalar(proba) && isvector(scale) && ~isscalar(scale)
+ proba = repmat(proba, length(scale), 1);
end
-if proba>1-2*eps()
- t = Inf;
+%
+% Compute inverse CDF
+%
+
+t = NaN(size(proba));
+t(proba<2*eps()) = 0;
+t(proba>1-2*eps()) = Inf;
+
+if all(~isnan(t))
return
end
-t = exp(log(scale)+log(-log(1-proba))/shape);
+i = (proba>=2*eps() | proba<=1-2*eps());
+
+t(i) = exp(log(scale(i))+log(-log(1-proba(i)))./shape(i));
return % --*-- Unit tests --*--
%@test:1
try
- x = wblinv(0, 1, 2);
- t(1) = true;
+ x = wblinv(0, 1, 2);
+ t(1) = true;
catch
t(1) = false;
end
@@ -76,8 +102,8 @@ T = all(t);
%@test:2
try
- x = wblinv(1, 1, 2);
- t(1) = true;
+ x = wblinv(1, 1, 2);
+ t(1) = true;
catch
t(1) = false;
end
@@ -94,16 +120,16 @@ shapes = [.1, 1, 2];
x = NaN(9,1);
try
- k = 0;
- for i=1:3
- for j=1:3
- k = k+1;
- x(k) = wblinv(.5, scales(i), shapes(j));
- end
- end
- t(1) = true;
+ k = 0;
+ for i=1:3
+ for j=1:3
+ k = k+1;
+ x(k) = wblinv(.5, scales(i), shapes(j));
+ end
+ end
+ t(1) = true;
catch
- t(1) = false;
+ t(1) = false;
end
if t(1)
@@ -126,44 +152,114 @@ x = NaN(9,1);
p = 1e-1;
try
- k = 0;
- for i=1:3
- for j=1:3
- k = k+1;
- x(k) = wblinv(p, scales(i), shapes(j));
- end
- end
- t(1) = true;
+ k = 0;
+ for i=1:3
+ for j=1:3
+ k = k+1;
+ x(k) = wblinv(p, scales(i), shapes(j));
+ end
+ end
+ t(1) = true;
catch
- t(1) = false;
+ t(1) = false;
end
if t(1)
- k = 1;
- for i=1:3
- for j=1:3
- k = k+1;
- shape = shapes(j);
- scale = scales(i);
- density = @(z) exp(lpdfgweibull(z,shape,scale));
- if debug
- [shape, scale, x(k-1)]
- end
- if isoctave
- s = quadv(density, 0, x(k-1),1e-10);
- else
- s = integral(density, 0, x(k-1));
- end
- if debug
- [s, abs(p-s)]
- end
- if isoctave
- t(k) = abs(p-s)<1e-9;
- else
- t(k) = abs(p-s)<1e-12;
- end
- end
- end
+ k = 1;
+ for i=1:3
+ for j=1:3
+ k = k+1;
+ shape = shapes(j);
+ scale = scales(i);
+ density = @(z) exp(lpdfgweibull(z,shape,scale));
+ if debug
+ [shape, scale, x(k-1)]
+ end
+ if isoctave
+ s = quadv(density, 0, x(k-1),1e-10);
+ else
+ s = integral(density, 0, x(k-1));
+ end
+ if debug
+ [s, abs(p-s)]
+ end
+ if isoctave
+ t(k) = abs(p-s)<1e-9;
+ else
+ t(k) = abs(p-s)<1e-12;
+ end
+ end
+ end
end
T = all(t);
%@eof:4
+
+%@test:5
+try
+ x = wblinv(.3, [1; .5], [2; 1.5]);
+ t(1) = true;
+catch
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isvector(x) && length(x)==2;
+end
+T = all(t);
+%@eof:5
+
+%@test:6
+try
+ x = wblinv(.3, [1; .5], 1.5);
+ t(1) = true;
+catch
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isvector(x) && length(x)==2;
+end
+T = all(t);
+%@eof:6
+
+%@test:7
+try
+ x = wblinv(.3, .5, [2; 1.5]);
+ t(1) = true;
+catch
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isvector(x) && length(x)==2;
+end
+T = all(t);
+%@eof:7
+
+%@test:8
+try
+ x = wblinv([.3;.1], .5, [2; 1.5]);
+ t(1) = true;
+catch
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isvector(x) && length(x)==2;
+end
+T = all(t);
+%@eof:8
+
+%@test:9
+try
+ x = wblinv([.3;.1], .5, 1.5);
+ t(1) = true;
+catch
+ t(1) = false;
+end
+
+if t(1)
+ t(2) = isvector(x) && length(x)==2;
+end
+T = all(t);
+%@eof:9