Issues with qbeta and pbeta I'm trying to get samples from truncated beta. I have written the following function 
based on inverse transform sampling method to do so, however, it seems that when CDF converges to 1 fast, pbeta and qbeta have some numerical issues.    
rTruncBeta <- function(n, alpha, beta, lbound, ubound) {
 F0 <- pbeta(lbound, alpha, beta)
 F1 <- pbeta(ubound, alpha, beta)
 U <- runif(n, F0, F1)
 qbeta(U, alpha, beta)
}

Suppose we need to generate samples from truncated beta between .39 and .43 with $\alpha = .1$ and $\beta = 57$.But sometimes this function generates 1 rather than generating between the two bounds. How to fix this issue (of course without sampling more and discarding the 1s)
> rTruncBeta(10,.1, 57,.39, .43)

[1] 0.3937340 0.3920652 0.3928667 1.0000000 0.4004606 0.3924584   
   1.0000000 0.3902923 0.3913201 0.4046171

 A: Note that pbeta(0.39,0.1,57) is about $1-3.6\times 10^{-15}$ and  pbeta(0.43,0.1,57) is about $1-7.0\times 10^{-17}$.
Near-total loss of numerical accuracy is the problem. 
Indeed, getting "1"s is far from your only issue here. If I generate a sample from your function and a better alternative in that region of the parameter space (see below)  and put into tb and tb2 respectively:
tbi <- tb[tb<.5]
qqplot(tb2[1:length(tbi)],tbi)

so as to compare the quantiles:

you can see that your function is far from continuous, taking only a few distinct values at the upper half of its interval.
You could put tests in your TruncBeta function to avoid such catastrophic loss of accuracy - if you check when you're in the far tail like that, you can always call something similar that retains a lot more significant figures (and so will look more like what you expect it to). For example, up in the extreme tail, try something like
rTruncBeta2 <- function(n, alpha, beta, lbound, ubound) {
 S0 <- pbeta(lbound, alpha, beta, lower.tail=FALSE)
 S1 <- pbeta(ubound, alpha, beta, lower.tail=FALSE)
 U <- runif(n, S1, S0)
 1-qbeta(U, beta, alpha)
}

