How to handle large numbers of output data sets from a simulation/sensitivity analysis?

Question

Somewhat related, but I think the question is distinct enough to justify a separate question.

As a bit of background, I come from a observational/statistical Epidemiology background, working with data that's been collected already, so even our large datasets tend to be in discrete, inseparable chunks. As such, I've never actually learned how to handle large amounts of data coming out of a simulation.

My problem is as follows:

I'm working on a compartmental model of an infection system, involving a system of ~10 differential equations and ~40 parameters shared among them. Most of these parameters aren't constants, but rather statistical estimates, and as such could be drawn from a distribution. That's one of the things I'd like to do - see how much the system varies purely as the result of uncertainty in the parameter estimates.

Which involves running the numerical solution to the model sampling from each distribution many many times to cover the parameter space. If it was just a stochastic simulation, I might be fine with outputting one huge, or thousands of small data files that I could churn through with a script. My current problem is how to manage the output data given I need to know what parameter values were drawn.

Right now, my very quick and dirty method is to output the parameter value alongside the numerical results, meaning that if I ran the system for 100 steps, I end up with a 100 row, 50 column data set - but where 40 of those columns are the same number repeated over and over again. That seems hugely wasteful, and is resulting in really large files.

Surely there's a better way to do this? Most of this is, currently, being implemented in Python.

Answer 1

You could try using something like a very basic Relational database. You could label every output file with a separate key, e.g. a sequential number, and then maintain a separate file which, in every row, contains the keys and the parameters used.

If you're processing the data automatically, you will have to use one level of indirection, but that still saves you quite a few table entries compared to your model.

On a side note, it sounds like you're doing parameter sensitivity analysis by just sampling the parameter space. If you're trying to minimize some function of the results of your ODE integration, you may want to look at approaches which also integrate the derivative of the system variables with respect to the parameters. I have recently published a paper on this (sorry for the shameless plug), which you might find interesting, if not for the references to other methods.

Answer 2

If you want to analyze the sensitivity of your outcome on changes in the input space, you would want something like Dakota or SUSA (which is described in this paper). These codes allow you to run your simulation as a black box a number of times while sampling the parameters from the probability distribution you assign to them. The output is a statistical result: you get a mean and standard deviation on your output value. SUSA will also rank the input parameters according to sensitivity etc.

In the case of ODEs, you might want to go for the Sundials package that has specific solvers that allow sensitivity analysis on parameters in the ODE. The Assimulo package has Python bindings for these solvers.

Answer 3

As the main goal is to compute afterwards some statistics, the correct approach is to decide on a genertous set of additive sufficient statistics for what you possibly want to compute afterwards, and then to update thes statistics with each batch. Then you never need to store old data, though you misgt store some as backup or for independent hypothesis testing.

The relevant statistics typically consist of a collection of numbers computed as the sum or max or min of selected functionals (key features, products of key features, key features in the presence of key indicators, etc.) together with the total number of items summed (maxxed, minned) over. This allows you to do at any instant correspondng means, covariance matrices, and probability tables of interest for subsequent statistical interpretation.

If you want to do some sort of 5-fold cross-validation, say, to validate your statistical conclusions, make 5 such collections, assigning each new item randomly to one of the five collections. Then combine the statistics for 4 out of the 5 collections to make predictions, and check their accuracy on the fifth one (each one taking this role in turn). This gives 5 accuracy estimates.