code: def _calculate_plotting(n, incidence):
          
          p = incidence / float(n)
          z = stats.norm.ppf(0.975)
          q = 1. - p
          ll = (
              (2 * n * p + 2 * z - 1) - z * np.sqrt(2 * z - (2 + 1 / n) + 4 * p * (n * q + 1))
          ) / (2 * (n + 2 * z))
          ul = (
              (2 * n * p + 2 * z + 1) + z * np.sqrt(2 * z + (2 + 1 / n) + 4 * p * (n * q - 1))
          ) / (2 * (n + 2 * z))
          return p, ll, ul
docs: Add confidence intervals to dichotomous datasets. From bmds231_manual.pdf, pg 124-5.
  
          LL = {(2np + z2 - 1) - z*sqrt[z2 - (2+1/n) + 4p(nq+1)]}/[2*(n+z2)]
          UL = {(2np + z2 + 1) + z*sqrt[z2 + (2-1/n) + 4p(nq-1)]}/[2*(n+z2)]
  
          - p = the observed proportion
          - n = the total number in the group in question
          - z = Z(1-alpha/2) is the inverse standard normal cumulative
                distribution function evaluated at 1-alpha/2
          - q = 1-p.
  
          The error bars shown in BMDS plots use alpha = 0.05 and so
          represent the 95% confidence intervals on the observed
          proportions (independent of model).