Suppose we have a classical statistical problem with canonical coordinates $\vec{q} = (q_1, q_2, \dots, q_n)$ and $\vec{p} = (p_1, p_2, \dots, p_n)$ such that they fulfill the usual Poisson brackets: \begin{align} \{ q_i, p_j \} & = \delta_{i,j} \\ \{ q_i, q_j \} & = 0 \\ \{ p_i, p_j \} & = 0 \\ \end{align} One can show that the time evolution of every dynamical quantity $f(\vec{q}, \vec{p}, t)$ is given by $$ \frac{d f}{d t} = \frac{\partial f}{\partial t} + \{ f, H \} = \frac{\partial f}{\partial t} + \sum_{i=1}^n \frac{\partial f}{\partial q_i} \dot{q_i} + \frac{\partial f}{\partial p_i} \dot{p_i} = \frac{\partial f}{\partial t} + (\nabla f) \cdot \vec{v} = \frac{\partial f}{\partial t} + \nabla (f \cdot \vec{v}) $$ with $\nabla = (\frac{\partial}{\partial q_1}, \dots, \frac{\partial}{\partial q_n}, \frac{\partial}{\partial p_1}, \dots, \frac{\partial}{\partial p_n})$, $\vec{v} = (\dot{q_1}, \dots, \dot{q_n}, \dot{p_1}, \dots, \dot{p_n})$ and $H$ the Hamiltonian of the system.
Liouville's theorem states that: $$ \int d^n p ~ ~ d^n q = \int d^n p' ~ ~ d^n q' $$ if $(\vec{q}~', \vec{p}~')$ and $(\vec{q}, \vec{p})$ are both canonical coordinates, e.g. related by a canonical transformation. So the phase space volume is a constant between systems which are described by canonical coordinates.
Now consider the phase space density $\varrho(\vec{q}, \vec{p}, t)$ which is the density of dynamically allowed trajectories at a given point $(\vec{q}, \vec{p})$ in phase space for a given instance of time $t$.
Liouville's equation reads: $$\frac{d \varrho}{d t} = 0$$ which (together with the equation for $f$) says that $\varrho$ is a (locally) conserved density in phase space. Because $\varrho \ge 0$ one can conclude that there are no sources of trajectories in phase space: trajectories do not start, end or cross.
Usually the Liouville equation is derived from Liouville's theorem. However I haven't seen such a derivation for which at some point it wasn't assumed that $\varrho$ is a (locally) conserved density. Hence, that reasoning is circular.
Do you know a non-circular derivation of Liouville's equation or is it indeed an axiom of classical statistical mechanics?