diff --git a/matlab/@dprior/admissible.m b/matlab/@dprior/admissible.m
index 0bc685dae68f734b0c724dddb35b863754b380b2..02c6d8300913f6eeea3617ee7a13e3b0d33de80c 100644
--- a/matlab/@dprior/admissible.m
+++ b/matlab/@dprior/admissible.m
@@ -1,4 +1,4 @@
-function b = admissible(o, d)
+function [b, lbid, ubid] = admissible(o, d)
 
 % Return true iff d is an admissible draw in a distribution characterized by o.
 %
@@ -8,9 +8,12 @@ function b = admissible(o, d)
 %
 % OUTPUTS
 % - b      [logical]  scalar.
+% - lbid   [integer]  vector, indices of the elements (parameters) violating the lower bound.
+% - ubid   [integer]  vector, indices of the elements (parameters) violating the upper bound.
 %
 % REMARKS
-% None.
+% - An admissible draw must lie within the lower and upper bounds of the truncated prior mass.
+% - if b==true, then lbid and ubid are empty vectors.
 %
 % EXAMPLE
 %
@@ -23,7 +26,7 @@ function b = admissible(o, d)
 %
 %   1
 
-% Copyright © 2023-2024 Dynare Team
+% Copyright © 2023-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -42,12 +45,23 @@ function b = admissible(o, d)
 
 b = false;
 
+if nargout>1
+    lbid = [];
+    ubid = [];
+end
+
 if ~isequal(length(d), length(o.lb))
     return
 end
 
 if all(d>=o.lb & d<=o.ub)
     b = true;
+    return
+end
+
+if nargout>1
+    lbid = find(d<o.lb);
+    ubid = find(d>o.ub);
 end
 
 return % --*-- Unit tests --*--
@@ -142,5 +156,101 @@ catch
 end
 
 T = all(t);
-
 %@eof:1
+
+%@test:2
+% 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);
+
+p3(9) = 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);
+    % Do simulations in a loop and estimate recursively the mean and the variance.
+    draw = o.deviate();
+    draw(1) = -.1;  % beta
+    draw(2) = -.1;  % gamma
+    draw(8) = 1.1;  % beta
+    draw(9) =  .1;  % (shifted gamma)
+    [b, lbid, ubid] = o.admissible(draw);
+    t(1) = ~b;
+    lbid
+    ubid
+catch
+    t(1) = false;
+end
+
+t(2) = isequal(lbid, [1; 2]);
+t(3) = isequal(ubid, 8);
+
+T = all(t);
+%@eof:2