```
In [1]:
```% ignore - internal setup
path('../scripts/', path)
set(0, 'DefaultLineLineWidth', 2)
set(0, 'DefaultAxesFontSize', 14);

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];

```
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

```
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');

```
```

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

```
```

```
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 [ ]:
```