function [rho,z] = f_G_update_VB(rho,linpred,residual,tau,ltau)

%% |------------------- Function to update cluster index G
%% Tractable Bayesian density regression via logit stick-breaking priors
%% Tommaso Rigon & Daniele Durante, Journal of Statistical Planning and Inference 211 (2021) 131–142

% Original Code at https://github.com/tommasorigon/LSBP/blob/master/src/functions.cpp
% This version is translated from R/C++. 
% Translated by Maximilian Schröder BI Oslo, Norges Bank.

% X1    = variables in regression
% X2    = variables responsible for mixing  
% gamma = parameters of log-odds of LSBP
% beta  = regression parameters in mixtures
% tau   = precision parameter

%% |---------- Initialisations
H      = size(residual,2); % number of clusters
n      = size(residual,1);    % number of observations

zl     = zeros(n,H);
z      = zeros(n,H);
cum_nu = zeros(H,1);   % cumulative product of stick-breaking weights

rho = [rho ones(n,1)];
eta        = 0;


    % cycle through observations
    for i = 1:n
        for h = 1:H-1
            cum_nu(1) = 1; 
            for l = 1:H
                % Positive contributions
                if l == h
                    if h == 1
                        zl(i, l) = 1;
                    else
                        indexes = 1:l-1; % 1, 2, ..., l-1
                        zl(i, l) = prod(1 - rho(i, indexes));
                    end
                elseif l > h
                    indexes = 1:l-1; % 1, 2, ..., l-1
                    indexes(indexes == h) = [];
                    if h == 1 && l == 2
                        zl(i, l) = -rho(i, l);
                    else
                        zl(i, l) = -rho(i, l) * prod(1 - rho(i, indexes));
                    end
                else
                    zl(i, l) = 0;
                end  
            end

            % Linear predictor
            eta = linpred(i, h) + sum(zl(i, :).* (ltau/2 - tau/2.* residual(i, :))); % check again
            rho(i, h) = 1 / (1 + exp(-eta));

            % Stick-breaking
            z(i, h) = rho(i, h) * cum_nu(h);
            cum_nu(h+1) = cum_nu(h) * (1 - rho(i, h));
        end
        z(i, H) = cum_nu(H);
    end
    rho(:, H) = [];
  

end

