I'm reading about statistical definition of entropy, which says
$$S=-k_B\sum_ip_i\ln p_i,\tag1$$
where $k_B$ is Boltzmann's constant, and $p_i$ is probability of $i$th state to be occupied. But in classical mechanics there's always an uncountable set of states with given energy, e.g. all states between time $t_0$ and $t_1$ for given initial conditions.
Shouldn't $(1)$ be changed to an integral instead of discrete sum? Why is it never (as it seems) done?
Let us start with an example, called Langevin paramagnetism, where the magnetic moment is described classically, as a vector in three dimensions. Calling $\vec\mu$ this moment, $\vec B$ the magnetic induction and $\theta$ the angle between $\vec\mu$ and $\vec B$. The probability density of the angle $\theta$ is $\rho(\theta)=\frac{1}{\mathcal Z}\exp(\beta\mu B\cos\theta)$, where $\mathcal Z$ is the partition function and can be seen as a normalization constant. Indeed we have (using $x=\cos\theta$) $$\mathcal Z=\int_0^\pi\rho(\theta)\sin\theta\,\mathrm d\theta =\int_{-1}^1\exp(\beta\mu B x)\mathrm dx=2\frac{\sinh(\beta\mu B)}{\beta\mu B}.$$ The average energy is $$E=-\frac{\partial \ln\cal Z}{\partial \beta}=-\mu B\left(\coth(\beta\mu B)-\frac1{\beta\mu B}\right)$$ The entropy is $$\begin{split}S&=-k_{\text B}\int_0^\pi \rho(\theta)\ln\rho(\theta)\,\sin\theta\,\mathrm d\theta=-k_{\text B}\frac1{\cal Z} \int_{-1}^1\exp(\beta \mu B x)\left(\beta\mu Bx-\ln\mathcal Z\right)\mathrm dx\\ &=k_{\text B}\ln\mathcal Z-k_{\text B}\beta\mu B\left(\coth(\beta\mu B)-\frac{1}{\beta \mu B}\right)=k_{\text B}\ln {\cal Z}+\frac{E}{T}. \end{split}$$ If one defines (as usual) the free energy $F$ as $F=-k_{\text B}T\ln\mathcal Z$, then the last expression is $F=E-TS$. This demonstrates that the entropy formula is valid using the expression $S=-k_{\text B}\int\rho\ln\rho$.
In this model, no underlying discretization is needed, but when one uses the definition of $S$ as an integral, the argument in the logarithm may have a physical unit, which gives a hint that the entropy is defined up to a constant. Moreover, discretizing the values of $\theta$ into segments of width $\pi/N$, we can define $p_n=\rho(n\frac\pi N)\frac\pi N$. The statistical entropy becomes $$S=-k_{\text B}\sum_n p_n\ln p_n=-k_{\text B}\sum_n \rho\left(n\frac\pi N\right)\frac \pi N\left[\ln\rho\left(n\frac\pi N\right)+\ln\frac\pi N\right]$$ This is a Riemann sum plus a number that goes to infinity when $N\to\infty$. The preceding formula becomes $$S=-k_{\text B}\int\rho\ln\rho\;\;-k_{\text B}\lim_{N\to\infty}\ln\frac\pi N.$$ As usually only the variations of entropy are used, the infinite constant is irrelevant and one may forget about it and use the entropy defined by the integral.
The main reason sum of countable set instead of integral over continuous volume is introduced in statistical physics in the context of classical physics is simplification of mathematics needed to make use of the probability theory to explain basics of information theory and derive the very notion of information entropy, given by the $p\ln p$ formula. Or use Lagrange multipliers to get Boltzmann distribution.
It is easier to think and reason about formula that contain sums, and even easier if sums contain finite number of terms, like for probability theory of a die. One can explain what information entropy and Boltzmann distribution is based on probability theory easily that way.
Once the principles of information theory/statistical physics are explained, one can reiterate with more involved integral formalism that is more in line with classical mechanics and phase space.
In quantum theory, old founders decided that probabilities should be ascribed to eigenvalues of an operator instead of to states; consequently, the formula with sums is often even if the states possible form an even higher cardinality set than phase space does.
