function [out] = VB_BDR_old_dc(y,X1,X2,H,init_method,prior_mix,prior_kern,maxiter,tol)
    

    %% |---------------- Preliminaries & Infer Dimensions
    n        = size(y,1);
    p_kernel = size(X1,2);
    p_mixing = size(X2,2);
    log2pi   = log(2*pi);   % pre-compute log of pi
    
    %% |---------------- Prior

    b_kernel  = zeros(1, p_kernel);
    B_kernel  = eye(p_kernel);
    P_kernel  = eye(p_kernel)/B_kernel; % inverse of B_mixing
    Pb_kernel = P_kernel*b_kernel';
    eig_Bkernel = eig(B_kernel);
    logdet_Bkernel = sum(log(eig_Bkernel));
    
    b_mixing  = zeros(1, p_mixing);
    B_mixing  = eye(p_mixing);
    P_mixing  = eye(p_mixing)/B_mixing; % inverse of B_mixing
    Pb_mixing = P_mixing*b_mixing';
    eig_Bmixing     = eig(B_mixing); 
    logdet_Bmixing = sum(log(eig_Bmixing));
    
    a_tau = 0.001;
    b_tau = 0.001;
    
    %% Prior Mixing
    % lambda_j ~ N(0,D)
    Ddiag_mix = 10*ones(p_mixing,H-1); 
    Dinv_mix = repmat(1*eye(p_mixing),1,1,H-1);

    if prior_mix == 1
    % machine learning prior in mixing 
        a_mix = 1.0e-4;
        b_mix = 1.0e-4; 
    elseif prior_mix == 2
        % Horseshoe prior for Phi
        Ab_mix   = (0.0001)^2;1;(0.0001)^2;
        tau_mix  = 0.1*ones(H,p_mixing);   1*ones(H,p_kernel);0.1*ones(H,p_mixing);      
        nu_mix   = 0.1*ones(H,p_mixing);   1*ones(H,p_kernel);0.1*ones(H,p_mixing);  
    else
        Dinv_mix   = repmat(1*eye(p_mixing),1,1,H-1);
    end
    

    
    %% Prior Kernel
    % lambda_j ~ N(0,D)
    Ddiag_kern = 10*ones(p_kernel,H); 
    Dinv_kern = repmat(1*eye(p_kernel),1,1,H);

    if prior_kern == 1
    % machine learning prior in mixing 
        a_kern = 1.0e-4;
        b_kern = 1.0e-4; 
    elseif prior_kern == 2
        % Horseshoe prior for Phi
        Ab_kern   = (0.0001)^2;1;(0.0001)^2;
        tau_kern  = 0.1*ones(H,p_mixing);   1*ones(H,p_kernel);0.1*ones(H,p_kernel);      
        nu_kern   = 0.1*ones(H,p_mixing);   1*ones(H,p_kernel);0.1*ones(H,p_kernel);  
    else
        Dinv_kern = repmat(1*eye(p_kernel),1,1,H);
    end


    %% |---------------- Initialization

    % initialize output structures
%     beta_mixing_draws = zeros(H-1,p_mixing,draws/nthin);
%     beta_kernel_draws = zeros(H,p_kernel,draws/nthin);
%     tau_draws         = zeros(H,draws/nthin);
%     n_cluster         = zeros(1,draws/nthin);
%     G_draw            = zeros(n,draws/nthin);
%     logpost           = zeros(1,draws+burnin);
    
    tau            = ones(1,H)./diff(quantile(y,[0.25,0.75]));
    ltau           = zeros(1,H); 
    a_tilde        = zeros(1,H); 
    b_tilde        = zeros(1,H);
    xi             = ones(n,H-1);
    residual       = zeros(n,H);
    linpred_mixing = zeros(n,H-1);
    linpred_kernel = zeros(n,H);
    rho            = zeros(n,H-1);
    mu_mixing      = zeros(H - 1, p_mixing);
    Sigma_mixing   = zeros(H - 1, p_mixing, p_mixing);
    mu_kernel      = zeros(H, p_kernel);
    Sigma_kernel   = zeros(H, p_kernel, p_kernel); 
    

    % intialize sampler
    if strcmp(init_method,'random')
        % Random initialization
        mu_mixing = -0.1 + (0.1 - -0.1) .* rand((H - 1), p_mixing);
        Sigma_mixing = zeros(H - 1, p_mixing, p_mixing);

        % Random and uniform starting probabilities
        rho = rand(n, H - 1);
        z = zeros(n, H);
        for i = 1:n
          z(i, :) = f_sb(rho(i, :));
        end
    elseif strcmp(init_method,'cluster')
        for hh = 2:H
            G   = kmedoids(y,hh,'Distance', 'euclidean');
            s   = silhouette(y, G);
            silhouetteAvg(hh-1) = mean(s);
        end
        [~, bestH] = max(silhouetteAvg);
        G   = kmedoids(y,bestH+1);
        z = f_model_matrix(ones(n,1),G);
        z = z+1e-6;
        z = z./sum(z,2);
        z = [z zeros(size(z,1),H-(bestH+1))+1e-6];
        for h = 1:H - 1
            rho(:, h) = z(:, h) ./ sum(z(:, h:H), 2);
        end
    elseif strcmp(init_method,'standard')
        if H==2
            rho = 0.99;
        else
            rho = linspace(0.99,0.01,H-1); 
        end
        z = repmat(f_sb(rho),n,1);
        rho = repmat(rho,n,1);
    elseif strcmp(init_method,'flat')
        rho  = 1/(H).*ones(1,H-1);
        z = repmat(f_sb(rho),n,1);
        rho = repmat(rho,n,1);  
    end
    
    % initialize pieces of the lower bound
    lower1 = zeros(H, 1);
    lower4 = zeros(H, 1);
    lower6 = zeros(H, 1);
    lower2 = zeros(H - 1, 1);
    lower3 = zeros(H - 1, 1);
    lower5 = zeros(H - 1, 1);
    lower7 = zeros(H - 1, 1);
    lowerbound = -Inf;
    
    L = zeros(maxiter,1);

    for r = 1:maxiter
        for h = 1:H
           % logistic regressions
           if h < H
               omega = tanh(xi(:, h) / 2) ./ (2 * xi(:, h));
               omega(isnan(omega)) = 1 / 4;
               
               % Parameters for beta_mixing
               Sigma_mixing(h, :, :) = inv((X2.*sqrt(omega))' * (X2.* sqrt(omega)) + Dinv_mix(:,:,h));  %inv((X2.*sqrt(omega))' * (X2.* sqrt(omega)) + P_mixing);
               mu_mixing(h, :) = squeeze(Sigma_mixing(h, :, :)) * (X2' * (rho(:, h) - 1/2) + Pb_mixing);
               linpred_mixing(:, h) = X2 * mu_mixing(h, :)';
               
               % Update prior
               if prior_mix == 1
                   a = a_mix + 0.5;
                   b = b_mix + 0.5*(mu_mixing(h,:).^2 + diag(squeeze(Sigma_mixing(h,:,:)))');
                   Ddiag   = b./a;
                   Dinv_mix(:,:,h) = diag(1./(Ddiag + 1e-10));
               elseif prior_mix == 2
                   d4              = sum(1./tau_mix(h,:)) + 1;
                   tau_mix(h,:)    = 1./(Ab_mix.*nu_mix(h,:)) + 0.5*(length(mu_mixing(h, :))+1)./d4;           
                   nu_mix(h,:)     = 1./Ddiag_mix(:,h) + 1./(Ab_mix.*tau_mix(h,:))';
                   Ddiag_mix(:,h)  = min(1e+2,0.5*(mu_mixing(h, :).^2 + diag(squeeze(Sigma_mixing(h,:,:)))') + 1./nu_mix(h,:));
                   Dinv_mix(:,:,h) = diag(1./(Ddiag_mix(:,h)+ 1e-10));
               end
               

               % Variational step
               xi(:, h) = sqrt(linpred_mixing(:, h).^2 + sum(X2 * squeeze(Sigma_mixing(h, :, :)) .* X2, 2));
           end
           
           % Parameters for posterior of mu and tau
           Sigma_kernel(h, :, :) = inv(tau(h) * (X1.*sqrt(z(:, h)))'*(X1.*sqrt(z(:, h))) + Dinv_kern(:,:,h)); %eye(p_kernel)/(tau(h) * (X1.*sqrt(z(:, h)))'*(X1.*sqrt(z(:, h))) + Dinv_kern(:,:,h));   % inv(tau(h) * (X1.*sqrt(z(:, h)))'*(X1.*sqrt(z(:, h))) + Dinv_kern(:,:,h));
           mu_kernel(h, :) = squeeze(Sigma_kernel(h, :, :)) * ((X1.* z(:, h).* tau(h))'* y + Pb_kernel);
        
            % Update prior
            if prior_kern == 1
               a = a_kern + 0.5;
               b = b_kern + 0.5*(mu_kernel(h,:).^2 + diag(squeeze(Sigma_kernel(h,:,:)))');
               Ddiag   = b./a;
               Dinv_kern(:,:,h) = diag(1./(Ddiag + 1e-10));
            elseif prior_kern == 2
               d4              = sum(1./tau_kern(h,:)) + 1;
               tau_kern(h,:)    = 1./(Ab_kern.*nu_kern(h,:)) + 0.5*(length(mu_kernel(h, :))+1)./d4;           
               nu_kern(h,:)     = 1./Ddiag_kern(:,h) + 1./(Ab_kern.*tau_kern(h,:))';
               Ddiag_kern(:,h)  = min(1e+2,0.5*(mu_kernel(h, :).^2 + diag(squeeze(Sigma_kernel(h,:,:)))') + 1./nu_kern(h,:));
               Dinv_kern(:,:,h) = diag(1./(Ddiag_kern(:,h)+ 1e-10));
            end
        
           linpred_kernel(:, h) = X1 * mu_kernel(h, :)';
           residual(:, h) = y.^2 - 2 * y .* linpred_kernel(:, h) + linpred_kernel(:, h).^2 + sum(X1 * squeeze(Sigma_kernel(h, :, :)) .* X1, 2);
           
           a_tilde(h) = a_tau + sum(z(:, h)) / 2;
           b_tilde(h) = b_tau + sum(z(:, h) .* residual(:, h)) / 2;
           tau(h) = a_tilde(h) / b_tilde(h);
           ltau(h) = psi(a_tilde(h)) - log(b_tilde(h));
        
        end
        
        %Mixture probabilities
        [rho,z] = f_G_update_VB(rho,linpred_mixing,residual,tau,ltau);
        rho = min(max(rho, 1e-16), 1 - 1e-16); % Added for numerical stability
       
        for h = 1:H
            lower1(h) = sum(z(:, h) .* (ltau(h)/2 - log2(pi)/2 - tau(h)/2 .* residual(:, h)));
            lower4(h) = -p_kernel * log2(pi)/2 - logdet_Bkernel/2 - 0.5 * (sum(diag(P_kernel * squeeze(Sigma_kernel(h, :, :)))) +...
                (mu_kernel(h, :) - b_kernel) * P_kernel * (mu_kernel(h, :) - b_kernel)') + ...
                a_tau * log(b_tau) - gammaln(a_tau) - b_tau * tau(h) + (a_tau - 1) * ltau(h);
            lower6(h) = -p_kernel * log2(pi)/2 - f_ldet(squeeze(Sigma_kernel(h, :, :)))/2 - ...
                 p_kernel/2 + a_tilde(h) * log(b_tilde(h)) - gammaln(a_tilde(h)) - b_tilde(h) * tau(h) + ...
                 (a_tilde(h) - 1) * ltau(h);
             
            if h < H
                lower2(h) = sum(rho(:, h) .* linpred_mixing(:, h) - (linpred_mixing(:, h) + xi(:, h)) / 2 + log(normcdf(xi(:, h))));
                lower3(h) = -p_mixing * log2(pi)/2 - logdet_Bmixing/2 - 0.5 * (sum(diag(P_mixing * squeeze(Sigma_mixing(h, :, :)))) + ...
                    (mu_mixing(h, :) - b_mixing) * P_mixing * (mu_mixing(h, :) - b_mixing)');

                lower5(h) = -p_mixing * log2(pi)/2 - f_ldet(squeeze(Sigma_mixing(h, :, :)))/2 - p_mixing/2;

                lower7(h) = sum(rho(:, h) .* log(rho(:, h)) + (1 - rho(:, h)) .* log(1 - rho(:, h)));
            end
        end
        
        lowerbound_new = sum(lower1) + sum(lower2) + sum(lower3) + sum(lower4) - sum(lower5) - sum(lower6) - sum(lower7);
        
        % Break the loop at convergence
        if abs(lowerbound_new - lowerbound) < tol
            disp(['Convergence reached after ', num2str(r), ' iterations.']);          
            break;
        end
        
        lowerbound = lowerbound_new;
        L(r)       = lowerbound;
        % Display status

        if mod(r, 50) == 0
            disp(['Lower-bound: ', num2str(round(lowerbound, 4)), ', iteration: ', num2str(r)]);
        end
        
        if r == maxiter
            warning(['Convergence has not been reached after ', num2str(r), ' iterations.']);
        end

    end

    % Output
    cluster = zeros(size(z, 1), 1);
    for i = 1:size(z, 1)
        [~, idx] = max(z(i, :));
        cluster(i) = idx;
    end
    
    out.beta_mixing  = mu_mixing;
    out.beta_kernel  = mu_kernel;
    out.tau          = tau;
    out.z            = z;
    out.cluster      = cluster;
    out.l            = L(1:r-1);
    out.r            = r;
    out.a_tilde      = a_tilde;
    out.b_tilde      = b_tilde;
    out.Sigma_mixing = Sigma_mixing;
    out.Sigma_kernel = Sigma_kernel;
    out.mu_mixing    = mu_mixing;
    out.mu_kernel    = mu_kernel;

end
