Assume that $\vec{z}(t)$, the state at time $t$ of a particle in a two-dimensional space, can be fully described by its position and velocity:
$\vec{z}(t) = [r_x(t)\ r_y(t)\ v_x(t)\ v_y(t)]$.
Furthermore, assume that given the initial conditions $\vec{z}(t_I)$ one can numerically integrate a system of ordinary differential equations to obtain the state at an arbitrary final time, $t_F$ or stop integration if an event is detected before reaching $t_F$. In that case the numerical integration process will report the event kind and event time.
A simple example would be the process of dropping a small ball over a plate with holes; the plate is placed above and parallel the ground: some initial conditions would lead to a trajectory that passes through one of the holes and reaches the ground. Other initial conditions would lead to trajectories that hit the plate before reaching the ground.
Given a nominal initial state $\vec{z}(t_I)$ and covariance matrix $\Sigma$, how would you make a statistical analysis of the possible outcomes?
The parameters of the statistical analysis could be, for example, the fraction of balls that reach the ground, or the mean velocity of the particles hitting the plate, or which of the holes lets the most balls pass through, and so on.
A simple approach generates random initial states by using the Cholesky decomposition of $\Sigma = LL^*$ together with a vector of normal variates $\vec{u}$ to obtain random initial states $\vec{w}$:
$\vec{w}(t_I) = \vec{z}(t_I) + L\vec{u}$
The ordinary differential equations are numerically integrated for each of many hundreds or thousands of $\vec{w}(t_I)$, the outcomes stored, and the statistical analysis is performed on the collection of stored outcomes: brute force Monte Carlo, and all I know.
I keep reading about Bayesian Inference, and Variance Reduction Techniques, and MCMC, and Antithetic Variates, and---for the life of me---I cannot apply them to this simple problem.
Take antithetic variates: I made an experiment where I used both
$\vec{w}(t_I) = \vec{z}(t_I) + L\vec{u}$
and
$\vec{w}(t_I) = \vec{z}(t_I) + L(-\vec{u})$
in the simulations. The results were exactly the same as in the Brute Force approach: same variance of the results; same confidence intervals (by the way, I compute the confidence intervals using the bootstrap).
I keep reading about Metropolis-Hastings, and Gibbs Sampling, and Importance Sampling, but I just cannot extrapolate their typical let's integrate a function! or let's approximate $\pi$! examples to my domain.
Could you point me out to techniques that would allow me to reduce the variance (as compared to the brute-force Monte Carlo), or a place where MCMC or Bayesian inference is used to estimate parameters associated with a dynamical system?
