Monte Carlo Simulation Example
The following example code is written in C#, but it's almost identical in C++, and very similar in F#, Visual Basic, VB.NET, Java, and MATLAB. Complete example source codein all these languages is included with the Solver Platform SDK.
// Example program calling the Monte Carlo Simulation engine.
// Given two uncertain variables X and Y with distributions:
// X = [automatically fitted to sample data] and Y = PsiNormal(1,1)
// Perform 1000 trials and calculate the functions: X + Y and X * Yprivate Engine_Action SimEvaluator_OnEvaluate(Evaluator evaluator)
{
// First, we ask for the current uncertain variable values
double [] Var = evaluator.Problem.VarUncertain.Value.Array;
m_Log.Dump("Evaluation x1 = " + Var[0] + ", x2 = " + Var[1] + ", Trial = " +
evaluator.Problem.Engine.Stat.FunctionEvals);
// We compute the function values given the current variable values
evaluator.Problem.FcnUncertain.Value[0] = Var[0] + Var[1]; // X+Y
evaluator.Problem.FcnUncertain.Value[1] = Var[0] * Var[1]; // X*Y
return Engine_Action.Continue;
}private void Example18()
{
m_Log.Clear();
try
{
int nvars = 2, nfcns = 2;
double [] correl = new double [] { 1.0, 0.5, 0.5, 1.0 };
double [] sample = new double [] // for distribution fitting
{-3.44,0.69,-0.11,1.25,0.73,-1.66,1.57,-1.64,-0.51,-0.19,
1.24,-0.96,1.40,-0.01,-0.77,0.91,-0.72,0.85,-0.20,0.59,
1.65,1.40,0.19,-1.38,-0.45,-0.71,1.25,-1.18,-0.33,2.33,
0.73,-0.75,-0.63,-1.06,-0.71,0.33,-0.82,0.43,0.10,0.88,
0.11,-0.88,0.02,-0.94,0.05,-0.06,-1.87,-2.03,-0.89,1.14};
using (Problem problem = new Problem(Solver_Type.Simulate, nvars, nfcns))
{
problem.Evaluators[Eval_Type.Sample].OnEvaluate += new
EvaluateEventHandler(SimEvaluator_OnEvaluate);
problem.Solver.NumTrials = 1000;
Distribution dist = new Distribution();
dist.InitSample(sample); // set up the sample
dist.AutoFit(true, 0.95); // continuous data, 95% confidence
if (dist.FitSuccess)
{
m_Log.Dump("Distribution Type = " + dist.Name);
for (int i=0; i< dist.NumParams;i++)
m_Log.Dump("Param[" + i + "] = " + dist.Params[i]);
}
else return; // fitting unsuccessful
problem.Model.VarDistribution[0] = dist; // assign the fitted distribution
problem.Model.VarDistribution[1] = // define an analytic distribution
new Distribution("PsiNormal(1,1)");
problem.Model.DistCorrelations = // define the correlations
new DoubleMatrix(Array_Order.ByRow, nvars, nvars, correl);
problem.Solver.Simulate(); // perform the simulation
m_Log.Dump("Status = " + problem.Solver.SimulateStatus);
// display simulation results
m_Log.Dump("\r\nStatistics on Uncertain Variables");
for (int i = 0; i < nvars; i++)
m_Log.Dump("\r\nVariable " + i);
m_Log.Dump("Minimum = " + problem.VarUncertain.Statistics.Minimum[i]);
m_Log.Dump("Maximum = " + problem.VarUncertain.Statistics.Maximum[i]);
m_Log.Dump("Mean = " + problem.VarUncertain.Statistics.Mean[i]);
m_Log.Dump("StdDev = " + problem.VarUncertain.Statistics.StdDev[i]);
m_Log.Dump("Variance = " + problem.VarUncertain.Statistics.Variance[i]);
m_Log.Dump("Skewness = " + problem.VarUncertain.Statistics.Skewness[i]);
m_Log.Dump("Kurtosis = " + problem.VarUncertain.Statistics.Kurtosis[i]);
m_Log.Dump("10th Percentile = " + problem.VarUncertain.Percentiles[i, 9]);
m_Log.Dump("90th Percentile = " + problem.VarUncertain.Percentiles[i, 89]);
m_Log.Dump("\r\nStatistics on Uncertain Functions");
for (int i = 0; i < nfcns; i++)
m_Log.Dump("\r\nFunction " + i);
m_Log.Dump("Minimum = " + problem.FcnUncertain.Statistics.Minimum[i]);
m_Log.Dump("Maximum = " + problem.FcnUncertain.Statistics.Maximum[i]);
m_Log.Dump("Mean = " + problem.FcnUncertain.Statistics.Mean[i]);
m_Log.Dump("StdDev = " + problem.FcnUncertain.Statistics.StdDev[i]);
m_Log.Dump("Variance = " + problem.FcnUncertain.Statistics.Variance[i]);
m_Log.Dump("Skewness = " + problem.FcnUncertain.Statistics.Skewness[i]);
m_Log.Dump("Kurtosis = " + problem.FcnUncertain.Statistics.Kurtosis[i]);
m_Log.Dump("10th Percentile = " + problem.FcnUncertain.Percentiles[i,9]);
m_Log.Dump("90th Percentile = " + problem.FcnUncertain.Percentiles[i,89]);
}
}
catch (SolverException ex)
{
m_Log.Dump("Exception " + ex.ExceptionType);
m_Log.Dump("Exception " + ex.Message);
}
}