Yeah, I ended up going sampling way. I used Beta-Binomial model to estimate the success rate and then used the parameter estimate to generate increasingly larger synthetic datasets. From each such dataset, I estimated the parameter and especially the 95% HDI width.
For Beta-Binomial model, I wonder whether I can quickly compute the confidence interval? The Beta distribution is conjugate to Binomial, and the posterior parameters can be obtained w/o any sampling.