Monte-Carlo simulation is a sometimes elegant (and sometimes crude) method for simulating complex systems. Parameters that affect the system are selected from random distributions and the system response to these values is then calculated. Repeating this process many times produces often useful information about the system. The method is especially useful for examining non-linear systems that are difficult to write equations for.
Monte-Carlo simulation proves useful in the following areas:
- Discrete-Event Simulation
- Systems Science
- Bayesian Inference
This list is by no means exhaustive. Please feel free to comment if your favorite application of Monte-Carlo simulation is missing.
MCS-libre is a C++ library to facilitate Monte-Carlo simulation that I wrote about ten years ago. I’m writing about its use now so that I can reference it when I demonstrate SWIG in a future post. The library is freely available at http://sourceforge.net/projects/mcs-libre. Usage information, with a full example, is given below:
Usage
Sampling
A class called mcsLibre facilites sampling during simulation. Begin by including the necessary header file:
#include "mcsLibre.h"
Then declare an object of type mcsLibre:
mcsLibre samplesource;
This initializes MCS-libre using the default 100 random number seeds. Users can provide an alternate set of seeds if they wish.
One can then sample from normal and uniform distributions across multiple random number streams.
double i = samplesource.normal(mean, stdev, stream); double j = samplesource.uniform(min, max, stream);
The stream argument used above is an integer that specifies which random number stream to generate values from. It is imperative that statistically independent variables sample from different streams!
Collecting Simulation Statistics
A class called mcsLibreCollectStats assists with collection of statistical information during simulation. Begin by including the necessary header file:
#include “mcsLibreCollectStats.h”
Then declare an object of type mcsLibreCollectStats as follows:
mcsLibreCollectStats stats_var1;
To add random samples of a variable to the collection object, use the add_value function:
stats_var1.add_value(30.3);
Retrieving Statistical Information
The mcsLibreCollectStats class provides functions for extracting statistical information out of its objects:
double m1 = stats_var1.get_mean_value(); // mean double s1 = stats_var1.get_stdev_value(); // standard deviation double mn1 = stats_var1.get_min_value(); // minimum value double mx1 = stats_var1.get_max_value(); // maximum value double n1 = stats_var1.get_number_samples(); // number of samples double sum1 = stats_var1.get_sum_of_samples(); // sum of samples
Example
The rest of this post demonstrates the use of MCS-libre with an example simulation.
Simulating Reynold’s Number
Reynold’s number is a non-dimensional value used to describe flow around an object in a fluid dynamics system. It is the ratio of inertial forces to viscous forces, and is described by the equation
where V is the velocity of the object relative to the fluid, x is the distance the fluid travels across the object, and nu is the kinematic viscosity of the fluid. In this example we sample values for V, x, and nu from different distributions and then calculate the Reynold’s number from the sampled values. We do this thousands of times, storing each result, to gain information about how the Reynold’s number value responds to the sampled variable space.
Full Implementation of Example
/* Demonstration of MCS-libre toolkit for computing propagation of error in the Reynolds number. */
#include <iostream>
// include mcsLibre simulation library
#include "mcsLibre.h"
// include mcsLibre statistical probe class
#include "mcsLibreCollectStats.h"
main()
{
std::cout << "Performing Monte-Carlo simulation..." << std::endl;
// create object for simulation
mcsLibre reynolds;
// create probe for data collection
mcsLibreCollectStats reynolds_stats;
// define parameter types
double Re, Re_mean, Re_stdev; // Reynolds number of flow
double nu, nu_mean, nu_stdev; // kinematic viscosity of fluid
double v, v_mean, v_stdev; // velocity of flow
double x, x_mean, x_stdev; // distance considered
// set mean and standard deviations from measurements
nu_mean = 0.0000157; nu_stdev = 0.000002; // m^2 / s
v_mean = 2; v_stdev = 0.1; // m / s
x_mean = 12; x_stdev = 1; // m
// simulation loop
int samples = 15000; // number of simulations
for (int loopvar = 0; loopvar < samples; loopvar++)
{
// simulate nu from stream 1
nu = reynolds.normal(nu_mean, nu_stdev, 1);
// simulate v from stream 2
v = reynolds.normal(v_mean, v_stdev, 2);
// simulate x from stream 3
x = reynolds.normal(x_mean, x_stdev, 3);
// calculate a sample Reynolds number
Re = (v * x) / nu;
// send sample value to data collection probe
reynolds_stats.add_value(Re);
}
// retrieve mean value from simulation
Re_mean = reynolds_stats.get_mean_value();
// retrieve standard deviation value from simulation
Re_stdev = reynolds_stats.get_stdev_value();
std::cout << std::endl << "Simulation of Error Propagation " << "in Reynold's Number" << std::endl;
std::cout << "Re = (v * x) / nu" << std::endl << std::endl;
std::cout << "Let v = normal(" << v_mean << ", " << v_stdev << "), x = normal(" << x_mean << ", " << x_stdev << ")," << std::endl;
std::cout << "and nu = normal(" << nu_mean << ", " << nu_stdev << ")." << std::endl << std::endl;
std::cout << "Simulating " << samples << " Reynold's numbers yields:" << std::endl << std::endl;
std::cout << "Sample Mean = " << Re_mean << " Sample Standard Deviation = " << Re_stdev << std::endl;
std::cout << "Sample Max = " << reynolds_stats.get_max_value();
std::cout << " Sample Min = " << reynolds_stats.get_min_value() << std::endl;
std::cout << "Sum of Samples = " << reynolds_stats.get_sum_of_values() << std::endl;
std::cout << "Number of Samples = " << reynolds_stats.get_number_samples() << std::endl;
std::cout << std::endl;
}
Example Results
