This is mostly a question about math and less about physics. In essence, it is the difference between discrete count and a continuous count:
For simplification, let’s consider just the linear part of the numerator. The discrete count is
$\sum_{n=0}^\varepsilon n=\frac{\varepsilon(\varepsilon+1)}{2}=\frac{\varepsilon^2}{2}+\frac{\varepsilon}{2} \tag{1}$
Using Gauss summation formula. The smooth (continuous) count is
$\int_{0}^{\varepsilon} x \mathrm{d}x=\frac{\varepsilon^2}{2} \tag{2}$
The discrete count includes an extra linear element that is missing in the continuous count. This extra term comes from the excess “area” atop of the distribution function (check the image below).
Discrete counting includes an additional term of $x$. PS: This is just an example. To calculate the actual results one must use the corresponding distribution formulas.
A discrete sum includes an excess “error” that the differentiation process eliminates by making the size of the package ($\Delta x$) infinitesimally small. In the end a different result is achieved when comparing integration and summation. The important part here is that both, the sum and the integral converge or diverge together. Meaning that there is a sum associated to an integral even though they may be different. In the current case, the integral is
$$\frac{\int_0^\infty \varepsilon e^{-\beta\varepsilon}d\varepsilon}
{\int_0^\infty e^{-\beta\varepsilon} d\varepsilon}
=\frac{1}{\beta} \tag{3}.$$
And the sum is
$$\frac{\sum_{n=0}^\infty n e^{-\beta n}}
{\sum_{n=0}^\infty e^{-\beta n}}
=\frac{1}{e^{\beta}-1} \tag{4}.$$
So, the same setup has 2 valid mathematical solutions. However, experimentally we get a set of data points that are better fitted with $(4)$ and we are oriented to take that choice. An important physical consequence from this choice (which is obvious from the counting name) is that energy must be made up of discrete packages or quantum.