What does conservation of probability mean in Classical Mechanics and why is it true? In the context of the Liouville equation, regularly the conservation of probability is invoked. (Of course, the overall probability is always conserved but this is a truism and not what is meant here. See also the discussion here.)
But what concretely does this means in the context of Classical Mechanics and why is it true?
Let's say that for some reason only three initial states are possible $A$, $B$, and $C$ and we can describe our system using the probability distribution:
\begin{align}
\rho(t=0,A) &= 0.7 \notag \\
\rho(t=0,B) &= 0.2 \notag \\
\rho(t=0,C) &= 0.1 
\end{align}
The total probability to find our system in a region $R$ which contains $A$ and $B$ but not $C$ is therefore $P(R,t=0)=90$%.
Now, as time passes on ($t\to t=t_1$) $R$ gets dragged by the Hamiltonian flow and becomes $\tilde R$. Moreover, our phase space points $A$, $B$, $C$ also get dragged and become $\tilde A$, $\tilde B$, $\tilde C$. The statement regularly used in the derivation of the Liouville equation is
$$ P(\tilde R, t= t_1) = P(R,t=0)=90 \% \, . $$
Is this correct because the time-evolution of the phase space points $A$, $B$, $C$ and the time-evolution of the region $R$ are both described by Hamilton's equation? (Formulated simpler: because we move them around equally as time passes on?)
And secondly, what does it imply for our concrete probability distribution? Liouville's theorem tells us the phase space volume is constant. And this in combination with the conservation of probability tells us that $\frac{d \rho}{dt}=0$. But does this tells us that while the probabilities at our original locations can be wildly different ($\frac{\partial \rho }{\partial t}\neq 0$), e.g.
\begin{align}
\rho(t=t_1,A) &= 0.3 \notag \\
\rho(t=t_1,B) &= 0.5 \notag \\
\rho(t=t_1,C) &= 0.2  
\end{align}
we certainly have
\begin{align}
\rho(t=t_1,\tilde A) &= 0.7 \notag \\
\rho(t=t_1,\tilde B) &= 0.2 \notag \\
\rho(t=t_1,\tilde C) &= 0.1 \quad ?
\end{align}
 A: Let me preface my answer by saying that I hope not to get into a lengthy discussion! The two links you provided seem to cover quite a lot of ground already.
I prefer to think of this in the conventional way: the quantity of interest is a probability density (per unit volume in phase space). The ensemble may be visualized as a very dense cloud of points in phase space, which allows a more concrete, physical, picture of this idea of probability density. The equations of motion dictate how this cloud of points evolves in time. To me, this is more useful than trying to assign a discrete probability to discrete states in phase space.
Then the idea is fairly simple. In a small sub-volume of phase space you can count the number of state points, divide this number by the subvolume, and you have the density. When this evolves in time, you can track the way the sub-volume evolves: it might be a small hypercube initially, but it will distort and change shape with time. Nonetheless, if you put a containing surface around it (imagine just constructing some kind of envelope), the volume inside the surface will not change. This is the key property of the Hamiltonian flow. We also know that no trajectories will enter or leave this volume. So both the volume, and the number of state points, will remain constant, and hence the density in that region of space (which follows the flow) remains constant. The analogy is often made with the flow of an incompressible fluid in 3D space: this is a bit misleading (because a fluid is typically incompressible because of physical interactions between the particles, and in our case we are dealing with completely independent representative points in phase space) but it does show the similarity between the equations governing the flows.
