```python
%pylab inline
pylab.rcParams['figure.figsize'] = (16.0, 8.0)
```

    Populating the interactive namespace from numpy and matplotlib


# Summarizing the Monte Carlo output

The result of a Monte Carlo simulation is a set of samples from the probability distribution associated with the measurand
$$ \{\mathbf{y}^{(k)},\,k=1,\ldots,M\} $$

The aim is to derive characteristic information about the measurand from this set:

1) best estimate

2) uncertainty associated with the best estimate

3) intervals/regions with a prescribed coverage probability

## Univariate measurand

1) best estimate
\begin{align}
y =& mean \{ y^{(k)}, k=1,\ldots,M\} \\
=& \frac{1}{M} \sum_{k=1}^M y^{(k)}
\end{align}
2) uncertainty associated with the best estimate
\begin{align}
u_y =& std\{ y^{(k)}, k=1,\ldots,M\} \\
=& \frac{1}{M-1} \sum_{k=1}^M (y-y^{(k)})^2
\end{align}

### Exercise 4.1

Draw randomly from the normal distribution with mean $\mu=1.3$ and standard deviation $\sigma=0.4$ and calculate best estimate and its uncertainy using 100, 200, 300, ..., 100 000 draws. Compare estimate and uncertainty with mean and standard deviation of the original distribution for the different number of draws.


```python
from scipy.stats import norm

```

Intervals with a prescribed coverage probability can be calculated from the Monte Carlo outcome as follows

1) Sort the Monte Carlo outcome in ascending order 
``` python
    sort(Y)
```
2) For propability $P$ calculate the corresponding fraction of Monte Carlo draws $q=int(PM)$

3) Set lower bound of interval as $r=int(0.5(M-q))$ for a probabilistically symmetrical interval

4) Calculate upper bound of interval as $r+q$

### Exercise 4.2

Draw randomly from the distribution from Exercise 4.1 and calculate the 95% probabilistally symetrical coverage interval from 1000, 10000 and 100000 draws. Compare the result to the exact 95% interval.


```python
from scipy.stats import norm

P = 0.95 # sought probability of coverage interval


```

## Multivariate measurand

1) best estimate
\begin{align}
\mathbf{y} =& mean \{ \mathbf{y}^{(k)}, k=1,\ldots,M\} \\
=& \frac{1}{M} \sum_{k=1}^M \mathbf{y}^{(k)}
\end{align}
2) uncertainty associated with the best estimate
\begin{align}
U_\mathbf{y} =& cov\{ \mathbf{y}^{(k)}, k=1,\ldots,M\} \\
=& \frac{1}{M-1} \sum_{k=1}^M (\mathbf{y}-\mathbf{y}^{(k)})(\mathbf{y}-\mathbf{y}^{(k)})^T
\end{align}

### Exercise 4.3

Draw randomly from the normal distribution with mean 
$$\mathbf{\mu}=\left( \begin{array}{c}
0.4 \\ -1.5
\end{array}\right)
$$ 
and covariance
$$
\Sigma=\left(\begin{array}{cc}
0.09 & -0.2 \\ -0.2 & 1.44
\end{array}\right)
$$
and calculate best estimate and its uncertainy using 1000, 10000 and 100 000 draws. Compare estimate and uncertainty with mean and covariance of the original distribution for the different number of draws.


```python
from scipy.stats import multivariate_normal


```

Regions with a prescribed coverage probability can be calculated from the multivariate Monte Carlo outcome as follows

1) Calculate the Cholesky decomposition of the sample covariance matrix $U_{\mathbf{y}}=\mathbf{LL}^T$

2) Transform the Monte Carlo outcomes
$$ \mathbf{y}_{(k)} = \mathbf{L}^{-1}(\mathbf{y}^{(k)}-\mathbf{y})$$
and sort according to the distance measure
$$ d^2_{(k)} = \mathbf{y}_{(k)}^T\mathbf{y}_{(k)} $$

3) calculate $k_P$ such that a fraction $P$ of all Monte Carlo outcomes satisfies $d_{(k)}<k_P$

This defines the ellipsoidal region $(\mathbf{\eta}-\mathbf{y})^TU_{\mathbf{y}}^{-1}(\mathbf{\eta}-\mathbf{y})<k^2_P$

For a bivariate normal distribution, the factor for a 95% coverage ellipsoidal region is given as the 95% quantile of the $\chi^2$ distribution with 2 degrees of freedom.

### Exercise 4.4

Calculate 100 000 random draws from the distribution from Exercise 4.3 and calculate the 95% coverage region. Compare to the true 95% coverage region.


```python
def calculate_ellipse(mu, Sigma, kP):
    vals, vecs = linalg.eigh(Sigma)
    order = vals.argsort()[::-1]
    vals = vals[order]
    vecs = vecs[:,order]
    theta = degrees(np.arctan2(*vecs[:,0][::-1]))
    width, height = kP * sqrt(vals)
    return width, height, theta

```


```python
from scipy.stats import multivariate_normal, chi2
from matplotlib.patches import Ellipse

mu = array([0.4, -1.5])
Sigma = array([[0.09, -0.2],[-0.2, 1.44]])


```

054b2ddc72872be47ecaf091beeb0902c63a0889
ipynb
Jupyter Notebook
.ipynb_checkpoints/05 Summarizing the Monte Carlo output-checkpoint.ipynb
TanjaDorst/MonteCarloHandsOn
bb240f2f301a476cd61f0064b7b29e30b9c0ccc0
MIT
2018-04-16T16:25:30.000Z
2018-04-16T16:25:30.000Z
.ipynb_checkpoints/05 Summarizing the Monte Carlo output-checkpoint.ipynb
TanjaDorst/MonteCarloHandsOn
bb240f2f301a476cd61f0064b7b29e30b9c0ccc0
MIT
2020-09-09T07:08:09.000Z
2020-09-14T07:36:59.000Z
.ipynb_checkpoints/05 Summarizing the Monte Carlo output-checkpoint.ipynb
TanjaDorst/MonteCarloHandsOn
bb240f2f301a476cd61f0064b7b29e30b9c0ccc0
MIT
2018-04-12T13:51:04.000Z
2020-09-10T08:30:32.000Z
Qwen/Qwen-72B
1. YES
2. YES