function [CRPS] = f_pred_CPRS(y,X1,X2,beta,gamma,tau)

    % 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);
    
    % compute weights
    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
    
    CRPS = 0;
    for mm = 1:H
        %if round(tau(mm),4) == 0
        %    CRPS = CRPS;
        %else
            CRPS = CRPS + pii(mm)*A_(y,X1(1,:)*gamma(mm,:)',sqrt(1/tau(mm)));
        %end
        for nn = 1:H
            CRPS = CRPS -0.5*pii(mm)*pii(nn)*A_(X1(1,:)*gamma(mm,:)',X1(1,:)*gamma(nn,:)',sqrt((1/tau(mm))+1/tau(nn)));
        end
    end
    
end    
    

function [A] = A_(y,mu,sig)
    A = 2*sig*normpdf(y,mu,sig) + (y-mu)*(2*normcdf(y,mu,sig)-1);
end

