I'm sorry you are having this problem. I am not sure what is going wrong so I have passed it on to the developer.   For the (alpha, beta) parameters in your problem, the beta PDF is highly singular at  x=1. In particular, the PDF at x=1 is undefined.  The CDF is very, very, steep as x --> 1.  Here is the graph:   proc iml; mean = 0.975; sdev = 0.10; alpha = (mean**2)*(1-mean)/(sdev**2)-mean; beta = (alpha)*(1-mean)/mean; print alpha beta; x = do(0, 1, 0.0001); p = cdf("beta", x, alpha, beta); title "Graph of CDF Function"; call series(x, p); I suspect the reason that you can't compute the quantile for a large number such as 0.7 or 0.9 is that the x value for those probabilities is too close to 1. The root-finding algorithm to find the quantile is failing for these extreme parameter values. Let's use the CDF function to compute the probabilities for which the quantile is of the form (1 - epsilon) for epsilon = 10^{-k}, k = 10, 11, ..., 16:   g = 10##(-10 : -16); p = cdf("beta", 1 - g, alpha, beta); lbls = putn(g, "BEST5."); print p[L="1 - 10^{-k}" F=8.6 c=lbls]; In standard double-precision arithmetic, the smallest number that you can add/subtract from unity is about 2.2E-16 (this is called "machine epsilon"). You can see that after about p=0.7, the quantile function will not be able to compute x because it will be too close to 1.  In practice, this problem occurs for even smaller values of p because finding the quantile is a root-finding problem which typically involves trying to bracket the root (and thus work very close to 1).   Presumably, you want to find a way to work around this issue.   1. If you want to graph the quantile function, just reverse the horizontal and vertical axes:   title "Graph of Quantile Function"; call series(p, x);    2. If you need to compute quantiles for other reasons, you can exploit a symmetry in the Beta distribution which enables you to always evalute x values that are near 0 instead of x values near 1. The symmetry is  cdf("beta", x, alpha, beta) = 1 - cdf("beta", 1-x, beta, alpha)  In terms of the quantile function, the x that satisfies x = QUANTILE("beta", x, alpha, beta) also satisfies x = 1 - QUANTILE("beta", 1-p, beta, alpha) In this way you can always evaluate values of x <= 0.5. For example, if you are iterating over a sequence of probability values, you can branch according to whether the probability is less than or greater than 0.5, like this: p = do(0.01, 0.99, 0.01); x = j(1, ncol(p), .); do i = 1 to ncol(p); if p[i] <= 0.5 then x[i] = quantile("beta", p[i], alpha, beta); else x[i] = 1 - quantile("beta", 1-p[i], beta, alpha); end; call series(p, x); ... View more