

In [1]:

    
% ignore - internal setup
path('../scripts/', path)
set(0, 'DefaultLineLineWidth', 2)
set(0, 'DefaultAxesFontSize', 14);

Monte Carlo

This is the simple Monte Carlo example you worked on in class.

Consider the following objective and constraint

\begin{align*} f(x) &= x_1^2 + 2x_2^2 + 3x_3^2\\ c(x) &= x_1 + x_2 + x_3 \le 3.5 \end{align*}

At the point:

$$x = [1, 1, 1]$$

the standard deviation in $x$ is (normally distributed)

$$\sigma_x = [0.0, 0.06, 0.2]$$

Compute the following:

Output statistics for $f$ (mean, standard deviation, histogram)
Reliability of $c$



In [2]:

    
fun = @(x) x(1)^2 + 2*x(2)^2 + 3*x(3)^2;
con = @(x) x(1) + x(2) + x(3);

x = [1.0, 1.0, 1.0];
sigma = [0.00, 0.06, 0.2];

We will use randn, which gives us a random number k sampled from a normal distribution. It is sampled from a unit normal with zero mean and a standard deviation of 1 so to translate to an arbitrary mean and standard deviation the random value will be $$ x = \mu + k \sigma $$ If we sample enough times, we should get reliable statistics. Let's put our sampling in a function because we don't yet know how many samples we need.



In [ ]:

    
function [muf, stdf, reliability, f, c] = stats(x, sigma, fun, con, n)

    f = zeros(n, 1);
    c = zeros(n, 1);

    % sample many times
    for i = 1:n
        x1 = x(1) + randn(1)*sigma(1);
        x2 = x(2) + randn(1)*sigma(2);
        x3 = x(3) + randn(1)*sigma(3);

        f(i) = fun([x1, x2, x3]);
        c(i) = con([x1, x2, x3]);
    end

    % mean
    muf = mean(f);

    % standard deviation
    stdf = std(f);

    % reliability
    reliability = nnz(c <= 3.5) / n;

end

Let's evaluate this function for different values of n (number of samples) to see how long it takes to converge.



In [4]:

    
nvec = logspace(1, 6, 20);
muvec = zeros(20, 1);
stdvec = zeros(20, 1);
relvec = zeros(20, 1);

for i = 1:20
    [muvec(i), stdvec(i), relvec(i), ~, ~] = stats(x, sigma, fun, con, round(nvec(i)));
end
    
figure();
semilogx(nvec, muvec, '-o');
xlabel('mu');

figure();
semilogx(nvec, stdvec, '-o');
xlabel('sigma');

figure();
semilogx(nvec, relvec, '-o');
xlabel('reliability');

It looks like it takes about 100,000 ($10^5$) simulations for the statistics to converge. Let's rerun that case and check our the histogram and statistics.



In [5]:

    
n = 1e5;
[mu, sd, rel, f, c] = stats(x, sigma, fun, con, n);
mu
sd
rel
figure()
histogram(f);









    



mu =

    6.1240


sd =

    1.2349


rel =

    0.9915

Notice that the distribution skews to the right even though all the inputs were symmetric normal distributions. Because the function is nonlinear (quadratic in this case), symmetric inputs don't lead to symmetric outputs. In particular squaring amplifies the effect of positive variations.

Monte Carlo with LHS

In our discussion on surrogate-based optimization we will about Latin Hypercube Sampling (LHS). Let's apply LHS to Monte Carlo to see if it speeds up convergence.



In [ ]:

    
function [muf, stdf, reliability, f, c] = stats(x, sigma, fun, con, n)

    f = zeros(n, 1);
    c = zeros(n, 1);

    % generate LHS points before hand
    lhs = lhsdesign(n, 2);  % n samples for 2 dimensions
    % there is no variation in sigma1 for this case.
    rand2 = icdf('Normal', lhs(:, 1), x(2), sigma(2));
    rand3 = icdf('Normal', lhs(:, 2), x(3), sigma(3));

    % can also use lhsnorm directly, but the above is more general for other distributions

    % sample many times
    for i = 1:n
        x1 = x(1);
        x2 = rand2(i);
        x3 = rand3(i);

        f(i) = fun([x1, x2, x3]);
        c(i) = con([x1, x2, x3]);
    end

    % mean
    muf = mean(f);

    % standard deviation
    stdf = std(f);

    % reliability
    reliability = nnz(c <= 3.5) / n;

end



In [ ]:

    
% run for different sizes of n

muLHS = zeros(20, 1);
stdLHS = zeros(20, 1);
relLHS = zeros(20, 1);

for i = 1:20
    [muLHS(i), stdLHS(i), relLHS(i), ~, ~] = statsLHS(x, sigma, fun, con, round(nvec(i)));
end



In [21]:

    
figure();
semilogx(nvec, muvec, '-o');
hold on;
semilogx(nvec, muLHS, '-o');
xlabel('mu');
ylim([5.8, 6.3])

figure(); 
semilogx(nvec, stdvec, '-o');
hold on;
semilogx(nvec, stdLHS, '-o');
xlabel('sigma');
ylim([1, 2])

figure(); 
semilogx(nvec, relvec, '-o');
hold on;
semilogx(nvec, relLHS, '-o');
xlabel('reliability');
ylim([0.99, 1])

Convergence is much faster!



In [ ]: