function [res] = f_LSBP_density_C_test(y,X1,X2,beta,gamma,tau)
    
    H = size(gamma,1);
    n = size(X1,1);
    
    % initialize variables
    pdf    = zeros(n,1);
    pii    = zeros(H,1);
    nu     = zeros(H,1);
    cum_nu = zeros(H,1);
    res    = zeros(1,length(y));
    
    for qp =1:length(y)
        % cycle through the observations
        for i = 1:n

            % first mixture component
            nu(1)     = 1/(1+exp(-(X2(i,:)*beta(1,:)')));
            pii(1)    = nu(1);
            cum_nu(1) = 1-nu(1);

            mean_    = X1(i,:)*gamma(1,:)';
            sd       = sqrt(1/(tau(1)));
            pdf(i)   = pii(1)*normpdf(y(qp),mean_,sd);

            % update other mixture components until H-1
            for h = 2:H-1
                nu(h)     = 1/(1+exp(-(X2(i,:)*beta(h,:)')));
                pii(h)    = nu(h) * cum_nu(h-1);
                cum_nu(h) = (1-nu(h))*cum_nu(h-1);

                mean_  = X1(i,:)*gamma(h,:)';
                sd     = sqrt(1/(tau(h)));
                pdf(i) = pdf(i) + pii(h)*normpdf(y(qp),mean_,sd);
            end

            % update the last mixture component
            nu(H)  = 1;
            pii(H) = nu(H)*cum_nu(H-1);

            mean_  = X1(i,:)*gamma(H,:)'; 
            sd     = sqrt(1/(tau(H)));
            pdf(i) = pdf(i) + pii(H)*normpdf(y(qp),mean_,sd);
        end
        res(qp) = mean(pdf);
    end
    
end

