If you try the commands
curve(beta(x, 0.01), from=169, to=172, n=1000)
curve(beta(0.01, x), from=169, to=172, n=1000)
you can see a discontinuity due to some Inf values returned by the beta function, maybe related to the Inf value of function gamma
value out of range in 'gammafn'
even though the Inf value returned by the gamma function is due to a to large result, the Inf value of the beta function is not justify.
The calculation was gammafn(a) * gammafn(b) / gammafn(a+b), and the first product overflowed, even though all the values are finite.
The order of calculation has been changed for R > 2.15.2.