function [beta_mixing_sim,beta_kernel_sim,tau_sim] = f_simul_VB(out_VB,R)
    % Posterior density estimate for the VB model
    p_mix     = size(out_VB.beta_mixing,2);
    p_kern    = size(out_VB.beta_kernel,2);
    H         = size(out_VB.beta_kernel,1); 
    beta_mixing_sim = zeros(R, H - 1, p_mix);
    beta_kernel_sim = zeros(R, H, p_kern);
    tau_sim = zeros(R, H);

    for h = 1:H
        if h < H
            [V,D] = eig(squeeze(out_VB.Sigma_mixing(h, :, :)), 'matrix');
            A1 = V * sqrt(D);
            beta_mixing_sim(:, h, :) = out_VB.mu_mixing(h, :) + (randn(R, p_mix) * A1);
        end
        [V,D] = eig(squeeze(out_VB.Sigma_kernel(h, :, :)), 'matrix');
        A1 = V * sqrt(D);
        if ~isreal(A1)
            [V,D] = eig((squeeze(out_VB.Sigma_kernel(h, :, :))+squeeze(out_VB.Sigma_kernel(h, :, :))')/2, 'matrix');
            A1 = V * sqrt(D);
        end
        beta_kernel_sim(:, h, :) = out_VB.mu_kernel(h, :) + (randn(R,p_kern) * A1);
        tau_sim(:, h) = gamrnd(out_VB.a_tilde(h), 1 ./ out_VB.b_tilde(h), [R, 1]);
    end
end

