# Monte Carlo Simulation of Complex Dynamical System

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?

OK, if I understand the problem correctly you have

• An unknown initial state $X_i$
• A deterministic, yet complicated function $f(X_i, \Theta)$ with parameters $\Theta$ that describes your dynamical system, i.e the transition from the initial state $X_i$ to the end state $X_e$

Question 1: How to find the statistical properties of the end state $X_e$, given some distribution of $X_i$? I am not aware of any general method that performs significantly better than simply drawing from $X_i$ and following the trajectories (what you call brute force). This is of course assuming it's impossible to analytically or numerically forward the probability distribution for $X_i$. I'd be interested to hear if such a method exists, but the normal solution to this problem in Bayesian statistics (e.g. to calculate the posterior predictive distribution) is simply forwarding the uncertainty by simulation as described above.

Question 2: How can I make inference about $X_i$ or $\Theta$? This one has a clear answer - if you cannot forward your end state analytically, it is nearly certain that you cannot calculate your likelihood function

$$p(X_e | X_i, \Theta )$$

directly. In this case, you have to resolve to a method that approximates the likelihood from simulation. The two main methods are

• Approximate Bayesian Computation (ABC)
• Synthetic Likelihood

Our review Hartig et al., 2011, Statistical inference for stochastic simulation models - theory and application should point you to the right literature.

Addition: If simulation time is a serious constraint, you may consider constructing an emulator of your dynamical system, to generate a fast, approximate solution of your ODEs. I'm not on top of the relevant literature though, so I would suggest that you ask a new question to find out what the most promising technique for your problem would be.

• Thanks for the info --- it does seem that I am stuck with simulating the system to obtain an empirical posterior distribution. – Escualo Feb 2 '16 at 17:44
• An approximate posterior, yes. But I would say that's good news: you can do everything you want, just at the cost of a bit of simulation time. 10 years ago you would have had no options at all. – Florian Hartig Feb 2 '16 at 20:36
• Indeed --- I've had the supercomputer running for about a week and it is still not done... that is decades of computing time! – Escualo Feb 3 '16 at 18:41
• If your problem that complicated to compute, you may consider using an emulator. I have expanded the question in that regard. Lots of additional assumptions though. – Florian Hartig Feb 3 '16 at 19:19

I'm almost sure I'm missing something, you talk about days of computer time for the bruteforce approach, while it seems a simple molecular-dynamics-like problem which you should be able to model with e.g. velocity-verlet in real time with millions of particles in a modern GPU.

For a statistical approach I would look into statistical mechanics, maybe you could model a microstate probability for your system, starting from e.g. a grand canonical ensemble with some additional constraint. Once you have a model probability you can calculate your means sampling from your distribution with Metropolis.

Maybe better suited for https://physics.stackexchange.com/ ?

• The system I am actually simulating is much more complex than the example I placed in the question :) – Escualo Feb 4 '16 at 13:59
• hehe knew I was missing something :) wanted to just add a comment but you it's not allowed to newcomers... Still I believe you could ask in physics about a statistical mechanics approach – pip Feb 4 '16 at 14:13