Computing the Radon-Nikodym Derivative by use of Sampling Trajectories

Question

Assume that you have two trajectories $y(t) \in \mathbb{R}^{n}$ and $z(t) \in \mathbb{R}^{n}$ ($t \in [0,\infty)$). Their time evolution is given by the ordinary differential equations of motion \begin{equation} \dot{y}(t) = f_y(y(t)) , \qquad y(0)=y_0 \end{equation} and \begin{equation} \dot{z}(t) = f_z(z(t)) , \qquad z(0)=z_0. \end{equation} Both are ergodic and sampling the normed probability density distributions \begin{equation} \rho_y(x ) = \lim_{T \to \infty} \frac1T \int_{0}^{T} \delta (y(t)-x) dt \geq 0 \forall x \in \mathbb{R}^{n} , \qquad \int_{\mathbb{R}^n} \rho_y(x) dx = 1 \end{equation} and \begin{equation} \rho_z(x ) = \lim_{T \to \infty} \frac1T \int_{0}^{T} \delta (z(t)-x) dt \geq 0 \forall x \in \mathbb{R}^{n} , \qquad \int_{\mathbb{R}^n} \rho_z(x) dx = 1 \end{equation} with the dirac-delta-function $\delta(x)$, so that $\int_{\mathbb{R}^n} \delta(x) = 1 $. Assume that $\rho_y(x)$ is absolutely continuous in $\rho_z(x)$ ($\rho_y(x) \prec\prec \rho_z(x)$). Radon-Nikodym says now, that the derivative $f$ exists, so that \begin{equation} \int_E d\rho_y = \int_E f d\rho_z \end{equation} holds true. By use of a function $h(x)$, i would like to compute now numerically and efficiently the quantity \begin{equation} \int_E h d\rho_y - \int_E h d\rho_z. \end{equation} Does someone have any practical tips for doing so? Also literature about how to cope in a practical questions with Radon-Nikodym derivatives would be kind.