Monte Carlo Experimental Design

Overview

Monte Carlo sampling is essentially a brute-force technique in which random samples are taken until confidence that the entire space has been sampled is satisfactory.

Random numbers are used to create sampling points in each direction.

Think of Monte Carlo ray-tracing: you send out a whole bunch of rays, each in random directions, and from the result you determine the radiative flux. Mathematically, you're performing an integration by randomly sampling the function you want to integrate, then adding up all of the random samples:

$ \int f(x) dx \approx \frac{1}{N} \sum_{i} f( x_i ) $

Explanation

Transforming Variables

For a distribution that is a function of $ m $ variables $ x_1, \dots, x_m $:

Each variable has its own range, $ \alpha_i \leq x_i \leq \beta_i $

This range must be converted to $ [0,1] $

$ \hat{x}_i = \frac{ x_i - \alpha_i }{ \beta_i - \alpha_i } $

so that $ x_i \in \left[ 0, 1 \right] \forall i = 1 \dots m $

Log Scale

If you're using a log scale, i.e. sampling logarithmically more at $ \alpha $ than $ \beta $

$ \hat{x}_i = \frac{ \log{(x_i)} - \log{(\alpha_i)} }{ \log{(\beta_i)} - \log{(\alpha_i)} } $

Histogram of samples in normal, non-transformed space. Samples are preferentially grouped toward lower end of range.	When the samples are plotted in the logarithm-transformed space, the preferential grouping disappears, as each bin now represents an order of magnitude, which drastically stretches the higher-end bins (so they include more samples) and squishes the lower-end bins (so they include fewer samples).

Selecting Samples

An $ m $-element vector of random numbers is generated to correspond with a single sample.

A number of sample points are selected, and the $ \hat{x}_i $ and corresponding $ x_i $ are stored

the function is evaluated at each input variable value $ x_i $

all random input vectors and corresponding output vectors are stored/saved

Running in Parallel

Because the code suited itself to being run in parallel, I was able to use the Matlab Parallel Computing toolbox by changing the for loop over all samples that evaluates the function at each sample.

help parfor

Example

Problem Information

For details about the problem, including the input uncertainty map, see Example Problem for Experimental Design

What I am trying to learn:

1. What does "true" function look like via MC sampling?
2. How many samples are required for MC?
3. How high a degree does the polynomial go to, and is that a function of the number of samples?
4. How does the polynomial compare to a Fourier analysis?

Plan: save the state after 10, 100, 1k, 10k samples and evaluate after each.

Code

Main article: Monte Carlo Experimental Design Matlab Code

Computing Response Surface

See Response Surface Methodology for general information on response surface methodology.

(Some specific resopnse surface results go here, for this particular case)

References

Matlab

• Statistics toolbox documentation: http://www.mathworks.com/help/toolbox/stats/bq_676m-2.html
  • Very helpful page on linear regression in general: goes through linear reg., multiple linear reg., multivariate linear reg., polynomial reg., multivariate polynomial reg. (response surfaces), etc.

• x2fx function documentation: http://www.mathworks.com/help/toolbox/stats/x2fx.html

• regstats function documentation: http://www.mathworks.com/help/toolbox/stats/regstats.html

• scatter3 (3D scatter plot) function documentation: http://www.mathworks.com/help/techdoc/ref/scatter3.html
  • Modifying axes properties: http://www.mathworks.com/help/techdoc/ref/axes_props.html
  • Creating axes graphics objects: http://www.mathworks.com/help/techdoc/ref/axes.html

• surf function documentation: http://www.mathworks.com/help/techdoc/ref/surf.html