function [var] = f_pred_var(X1, X2, beta, gamma, tau)

    % infer dimensions
    H = size(gamma,1);
    n = size(X1,1);
    
    % initialize variables
    pred   = zeros(n,1);
    var    = zeros(n,1);
    pii    = zeros(H,1);
    nu     = zeros(H,1);
    cum_nu = zeros(H,1);
    
    % 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);
        
        pred(i) = pii(1)* X1(i,:)*gamma(1,:)';
        
        % 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);

            pred(i) = pred(i) + pii(h)* X1(i,:)*gamma(h,:)';
        end
    
        % update the last mixture component
        nu(H)  = 1;
        pii(H) = nu(H)*cum_nu(H-1);

        pred(i) = pred(i) + pii(H)* X1(i,:)*gamma(H,:)';
        
        % variance part
        var(i) = pii(1)* (X1(i,:)*gamma(1,:)'-pred(i))^2 + 1/tau(1);
        
        for h = 2:H-1
            var(i) = var(i) + pii(h)* (X1(i,:)*gamma(h,:)'-pred(i))^2 + 1/tau(h);
        end
        
        var(i) = var(i) + pii(H)* (X1(i,:)*gamma(H,:)'-pred(i))^2 + 1/tau(H);
        
    end
    
end

