docs: def _nbinom_ztrunc_p(mu, k_agg):
          """ Calculates p parameter for truncated negative binomial
  
          Function given in Sampford 1955, equation 4
  
          Note that omega = 1 / 1 + p in Sampford
          """
  
          p_eq = lambda p, mu, k_agg: (k_agg * p) / (1 - (1 + p)**-k_agg) - mu
  
          # The upper bound needs to be large. p will increase with increasing mu
          # and decreasing k_agg
          p = optim.brentq(p_eq, 1e-10, 1e10, args=(mu, k_agg))
          return p