To follow up: have a look at pymc3.distributions.dist_math.incomplete_beta which calls power-series and continued-fraction expansions for the incomplete beta function (beta CDF). As these are power expansions, they use scan; and one or both of these has a broken gradient.
You might unwrap the scan to some term of the approximation that is good enough; and this should deal with the problem.
Alternatively you can implement the gradients directly rather than using automatic differentiation - as power series are easily differentiated. Continued fractions are nastier, but sum-product formulas exist for their derivatives.