Calculates p parameter for truncated negative binomial

        Function given in Sampford 1955, equation 4

        Note that omega = 1 / 1 + p in Sampford

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