function p = draw(o)

% Return a random draw from the prior distribution.
%
% INPUTS
% - o    [dprior]
%
% OUTPUTS
% - p    [double]   m×1 vector, random draw from the prior distribution (m is the number of estimated parameters).
%
% REMARKS
% None.
%
% EXAMPLE
%
% >> Prior = dprior();
% >> Prior.add('alpha', 'shape', 'beta', 'mean', 0.5, 'std', 0.1);
% >> d = Prior.draw()

% Copyright © 2023-2026 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/>.

MAXITER = 1e6;
p = NaN(rows(o.lb), 1);
% Dispatch table: {name, index_vector, draw_function}.
% Each draw function takes an index vector and returns random draws.
families = {};
if o.has('uniform')
    families{end+1} = {'uniform', o.iduniform, ...
                       @(id) rand(length(id),1).*(o.p4(id)-o.p3(id)) + o.p3(id)};
end
if o.has('gaussian')
    families{end+1} = {'gaussian', o.idgaussian, ...
                       @(id) randn(length(id),1).*o.p7(id) + o.p6(id)};
end
if o.has('gamma')
    families{end+1} = {'gamma', o.idgamma, ...
                       @(id) gamrnd(o.p6(id), o.p7(id)) + o.p3(id)};
end
if o.has('beta')
    families{end+1} = {'beta', o.idbeta, ...
                       @(id) (o.p4(id)-o.p3(id)).*betarnd(o.p6(id), o.p7(id)) + o.p3(id)};
end
if o.has('invgamma1')
    families{end+1} = {'inverse gamma type 1', o.idinvgamma1, ...
                       @(id) sqrt(1./gamrnd(o.p7(id)/2, 2./o.p6(id))) + o.p3(id)};
end
if o.has('invgamma2')
    families{end+1} = {'inverse gamma type 2', o.idinvgamma2, ...
                       @(id) 1./gamrnd(o.p7(id)/2, 2./o.p6(id)) + o.p3(id)};
end
if o.has('weibull')
    families{end+1} = {'weibull', o.idweibull, ...
                       @(id) wblrnd(o.p7(id), o.p6(id)) + o.p3(id)};
end
for k = 1:length(families)
    fname = families{k}{1};
    id = families{k}{2};
    drawfn = families{k}{3};
    p(id) = drawfn(id);
    oob = find((p(id) > o.ub(id)) | (p(id) < o.lb(id)));
    iter = 0;
    while ~isempty(oob)
        iter = iter + 1;
        if iter > MAXITER
            error('dprior::draw: rejection sampling failed after %d iterations for %s distribution (parameters: %s). The prior bounds may be too tight.', ...
                  MAXITER, fname, mat2str(id(oob)))
        end
        p(id(oob)) = drawfn(id(oob));
        oob = find((p(id) > o.ub(id)) | (p(id) < o.lb(id)));
    end
end

return % --*-- Unit tests --*--

%@test:1
% Build a 14-parameter prior (2 of each distribution family) with add(),
% then verify that sample mean/variance match the theoretical moments.
prior_trunc = 1e-10;
Prior = dprior();
for rep = 1:2
    Prior.add(sprintf('p%d_beta',      rep), 'shape', 'beta',       'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_gamma',     rep), 'shape', 'gamma',      'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_normal',    rep), 'shape', 'normal',     'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_invgamma1', rep), 'shape', 'inv_gamma',  'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_uniform',   rep), 'shape', 'uniform',    'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_invgamma2', rep), 'shape', 'inv_gamma2', 'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
    Prior.add(sprintf('p%d_weibull',   rep), 'shape', 'weibull',    'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
end

ndraws = 1e5;
m0 = Prior.p1;
v0 = diag(Prior.p2.^2);

try
    for i = 1:ndraws
        draw = Prior.draw();
        m1 = m0 + (draw-m0)/i;
        m2 = m1*m1';
        v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i;
        m0 = m1;
    end
    t(1) = true;
catch
    t(1) = false;
end

if t(1)
    t(2) = all(abs(m0-Prior.p1)<3e-3);
    t(3) = all(all(abs(v0-diag(Prior.p2.^2))<5e-3));
end
T = all(t);
%@eof:1
