% quantile
function q = quantile(pvalue, cdf, xord)
% q = quantile(pvalue, cdf, xord);
% pvalue ~ 1 x p gives prob values 
% cdf ~ n x h gives n cdf functions.  Must increase across rows.
% xord ~ 1 x h gives x-values for cdf.  Must be increasing.
% q(i,j) ~ n x p interpolates x-quantiles off of cdf(i,:) at pvalue(j)
%   q(i, 1) = xord(1) if pvalue(1) <= cdf(i,1)
%   q(i, p) = xord(p) if pvalue(p) >= cdf(i,h)

p = length(pvalue);
n = size(cdf, 1);
h = size(cdf, 2);
q = NaN(n, p);
for i = 1:n
    for j = 1:p
        k = cdf(i,:) < pvalue(j);
        k = sum(k');
        k1 = k+1;
        if isnan(cdf(i,1));
            q(i,j) = NaN;
        elseif k == 0;
            q(i,j) = xord(1);
        elseif k == h;
            q(i,j) = xord(h);
        else
            q(i,j) = xord(k) + (xord(k1)-xord(k))* ...
                (pvalue(j)-cdf(i,k))/(cdf(i,k1)-cdf(i,k));
        end
    end
end
% end of function quantile