code: def fit_richness(self, atol=1.e-3, maxiter=50):
          if np.isnan(self.u).any():
              logger.warning("NaN signal probability found")
              return 0., 0., None
          if not np.any(self.u):
              logger.warning("Signal probability is zero for all objects")
              return 0., 0., None
          if self.f == 0:
              logger.warning("Observable fraction is zero")
              return 0., 0., None
          richness = np.array([0., 1./self.f, 10./self.f])
          loglike = np.array([self.value(richness=r) for r in richness])
          found_maximum = False
          iteration = 0
          while not found_maximum:
              parabola = ugali.utils.parabola.Parabola(richness, 2.*loglike)
              if parabola.vertex_x < 0.:
                  found_maximum = True
              else:
                  richness = np.append(richness, parabola.vertex_x)
                  loglike  = np.append(loglike, self.value(richness=richness[-1]))
                  if np.fabs(loglike[-1] - np.max(loglike[0: -1])) < atol:
                      found_maximum = True
              iteration+=1
              if iteration > maxiter:
                  logger.warning("Maximum number of iterations reached")
                  break
          index = np.argmax(loglike)
          return loglike[index], richness[index], parabola