function G = r2tcdf(x, muYZ, SigmaYZ, df)
% returns G(x) = CDF of ratio of 2 Student t's X = Y/Z, with DF dof.  
% J. Huston McCulloch Feb. 2010
% x ~ N X 1  selected x values
% muYZ ~ 2 X 1
% SigmaYZ ~ 2 X 2 (same for all x)
% G ~ N x 1 cdf values
% set W = Zx-Y
muWZ = muYZ;
SigmaWZ = SigmaYZ;
n = size(x,1);
G = NaN(n,1);
for i = 1:n
    if (isnan(x(i)) || isnan(muYZ(1)*muYZ(2))) 
      x(i)=NaN;    % mvncdf deals with NaNs, but mvtcdf doesn't!
    else
    muWZ(1) = x(i) * muYZ(2) - muYZ(1);
    SigmaWZ(1,1) = (x(i)^2)*SigmaYZ(2,2) + SigmaYZ(1,1) - 2*x(i)*SigmaYZ(1,2);
    SigmaWZ(1,2) = x(i) * SigmaYZ(2,2) - SigmaYZ(1,2);
    SigmaWZ(2,1) = SigmaWZ(1,2);
    rho = SigmaWZ(2,1)/sqrt(SigmaWZ(1,1)*SigmaWZ(2,2));
    if abs(rho) >= 1
        rho
    end
    CWZ = toeplitz([1 rho]);
    zWZ = muWZ./sqrt(diag(SigmaWZ));
       % mvtcdf will rescale Sigma to Correlations, but doesn't rescale 
       % arguments! 
    G(i) = mvtcdf(-zWZ', CWZ, df) + mvtcdf(zWZ', CWZ, df);
    end
end
% end of r2tcdf