function [cdf] = f_pred_cdf(X1,X2,beta,gamma,tau,threshold)

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

            cdf(i) = cdf(i) + pii(h)*normcdf(threshold,X1(i,:)*gamma(h,:)',sqrt(1/tau(h)));
        end
    
        % update the last mixture component
        nu(H)  = 1;
        pii(H) = nu(H)*cum_nu(H-1);

        cdf(i) = cdf(i) + pii(H)*normcdf(threshold,X1(i,:)*gamma(H,:)',sqrt(1/tau(H)));
    end
end

