function new_prior = laplace_approximation(o, prior_trunc)

% Return a new @dprior object where every prior has been replaced by a
% Gaussian with the same mean (p1) and standard deviation (p2) as the
% original. Bounds on the prior support are dropped (p3 = -Inf, p4 = Inf);
% truncation (lb, ub) is recomputed for the given prior_trunc.
%
% This implements the Laplace prior approximation used for penalized
% method-of-moments estimation, where only the first two moments of the
% prior are retained.
%
% INPUTS
% - o            [dprior]  input prior object.
% - prior_trunc  [double]  scalar in [0, 1], probability mass to remove in
%                          each tail when recomputing (lb, ub). Defaults to 0.
%
% OUTPUTS
% - new_prior    [dprior]  Gaussian-Laplace approximation of o.

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

if nargin < 2
    prior_trunc = 0;
end

new_prior = dprior();

n = length(o.p1);
new_prior.p1   = o.p1;
new_prior.p2   = o.p2;
new_prior.p3   = -Inf(n, 1);
new_prior.p4   =  Inf(n, 1);
new_prior.p5   = o.p1;
new_prior.p6   = o.p1;
new_prior.p7   = o.p2;
new_prior.p11  = o.p11;
new_prior.name = o.name;

new_prior.idgaussian  = (1:n)';
new_prior.idbeta      = [];
new_prior.idgamma     = [];
new_prior.idinvgamma1 = [];
new_prior.iduniform   = [];
new_prior.idinvgamma2 = [];
new_prior.idweibull   = [];

bounds = prior_bounds(new_prior, prior_trunc);
new_prior.lb = bounds.lb;
new_prior.ub = bounds.ub;

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

%@test:1
% Verify that laplace_approximation collapses every distribution family to a
% Gaussian with the original mean/std, on a 7-parameter prior built with add().
prior_trunc = 1e-10;
Prior = dprior();
Prior.add('p_beta',      'shape', 'beta',       'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
Prior.add('p_gamma',     'shape', 'gamma',      'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
Prior.add('p_normal',    'shape', 'normal',     'mean', 0.0, 'std', 1.0, 'prior_trunc', prior_trunc);
Prior.add('p_invgamma1', 'shape', 'inv_gamma',  'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
Prior.add('p_uniform',   'shape', 'uniform',    'lower_bound', 0, 'upper_bound', 1, 'prior_trunc', prior_trunc);
Prior.add('p_invgamma2', 'shape', 'inv_gamma2', 'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);
Prior.add('p_weibull',   'shape', 'weibull',    'mean', 0.4, 'std', 0.2, 'prior_trunc', prior_trunc);

try
    Laplace = Prior.laplace_approximation(prior_trunc);
    t(1) = true;
catch
    t(1) = false;
end

if t(1)
    t(2) = length(Laplace) == 7;
    t(3) = isequal(Laplace.idgaussian, (1:7)');
    t(4) = isempty(Laplace.idbeta);
    t(5) = isempty(Laplace.idgamma);
    t(6) = isempty(Laplace.idinvgamma1);
    t(7) = isempty(Laplace.iduniform);
    t(8) = isempty(Laplace.idinvgamma2);
    t(9) = isempty(Laplace.idweibull);
    % Means and stddevs preserved.
    t(10) = isequal(Laplace.p1, Prior.p1);
    t(11) = isequal(Laplace.p2, Prior.p2);
    % Hyperparameters of a Normal are (mean, std).
    t(12) = isequal(Laplace.p6, Prior.p1);
    t(13) = isequal(Laplace.p7, Prior.p2);
    % Mode equals mean for a Normal.
    t(14) = isequal(Laplace.p5, Prior.p1);
    % Support is unbounded.
    t(15) = all(Laplace.p3 == -Inf);
    t(16) = all(Laplace.p4 == Inf);
    % Names preserved.
    t(17) = isequal(Laplace.name, Prior.name);
end
T = all(t);
%@eof:1
