docs: def fit_mle(self, data, b=None):
          """%(super)s
          In addition to data, requires ``b``, the upper limit of the
          distribution.
          """
          # Take mean of data as MLE of distribution mean, then calculate p
          mu = np.mean(data)
          if not b:
              b = np.sum(data)
          p = _geom_solve_p_from_mu_vect(mu, b)
  
          # Just return float, not len 1 array
          if len(np.atleast_1d(p)) == 1:
              return float(p), b
          else:
              return p, b