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