Add method to compute or modify (in place) prior bounds.

matlab/@dprior/bounds.m

+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

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

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

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