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

    Populating the interactive namespace from numpy and matplotlib


# Drawing from multivariate distributions

Draws from a multivariate normal distribution
$$ N(\mathbf{\mu}, \mathbf{\Sigma}) $$
can be generated by 

1) Calculate the Cholesky decomposition $\mathbf{\Sigma} = \mathbf{R}^T\mathbf{R}$
``` python
    from numpy.linalg import cholesky
    R = cholesky(Sigma)
```

2) Generate standard normally distributed values $\mathbf{Z}$

3) Evaluate $$ \mathbf{X} = \mathbf{\mu 1} + \mathbf{RZ}$$


### Exercise 2.1

Draw 1000 samples from the bivariate distribution
$$ N\left(\left( \begin{array}{c}
0.2 \\ -1.0
\end{array}\right),
\left(\begin{array}{cc}
0.01 & -0.014 \\ -0.014 & 0.04
\end{array}\right)
\right)
$$


```python
# solution for exercise 2.1

draws = 1000

# mean
mu = array([0.2, -1.0])
# covariance matrix
Sigma = array([[0.01, -0.014],[-0.014, 0.04]])
# Cholesky factor
R = linalg.cholesky(Sigma)
# standard normal distributed values (dimension 2x1000)
Z = random.randn(2,draws)
# transformation to non-standard normal
X = dot(mu[:,newaxis], ones((1, draws))) + dot(R, Z)
# transpose to match output of built-in method
X = X.T
```


```python
figure()
scatter(X[:,0], X[:,1], 100, edgecolor="none");
```

Draws from the multivariate normal distribution can more easily generated using the built-in **scipy** functions

```python 
    from scipy.stats import multivariate_normal
    X = multivariate_normal(mu, Sigma, size=1000)
```
**Note**
Scipy *rvs* functions return arrays of shape (number of draws, size of mean)

### Exercise 2.2

Repeat Exercise 2.1 with the built-in scipy function and compare the results.


```python
# solution for exercise 2.2
from scipy.stats import multivariate_normal

Xb = multivariate_normal.rvs(mu, Sigma, size=draws)

figure()
scatter(Xb[:,0],Xb[:,1],100,color="b", edgecolor="none",label="built-in function")
scatter(X[:,0],X[:,1],100, color="r", edgecolor="none",label="formula",alpha=0.2)
legend();
```

## Draws using a copula function

In many practical cases knowledge about the input quantities is available in terms of their individual distributions and a correlation coefficient. This is insufficient to assign a unique multivariate distribution. Therefore, a copula function can be defined
$$
C(\mu_1,\ldots,\mu_N) = \mathbb{P} \left[ X_1\leq G_{X_1}^{-1}(\mu_1)\ldots,X_N\leq G_{X_N}^{-1}(\mu_N) \right]
$$

### Example copula functions

* all input quantities are mutually independent
$$ C(\mu_1,\ldots,\mu_N) = \prod_{k=1}^N \mu_k $$

* the input quantities are correlated with $\rho\equiv 1$
$$ C(\mu_1,\ldots,\mu_N) = \min_{k} \mu_k $$

* two input quantities are correlated with $\rho$
$$ C(\mu_1,\mu_2) = F_2(G_{X_1}^{-1}(\mu_1),G_{X_2}^{-1}(\mu_2),\rho) $$

The copula can be used to incorporate the correlation coefficient and the individual distributions $g_{X_i}$ to formally define a multivariate distribution. 

#### Example
Input quantities $X_1,X_2$ with correlation coefficient $\rho$ and
\begin{align}
X_1 \sim & N(\mu, \sigma) \\
X_2 \sim & U(a, b)
\end{align}

Use bivariate normal copula function:

1) Draw from bivariate standard normal distribution
$$ z \sim N\left(\mathbf{0}, \left(\begin{array}{cc}
1.0 & \rho \\ \rho & 1.0
\end{array}\right) \right)
$$





2) Evaluate cumulative distribution function of the copula
\begin{align}
\zeta_1 =& G_N(z_1) \\
\zeta_2 =& G_N(z_2) 
\end{align}

3) Evaluate inverse cumulative distribution functions
\begin{align}
x_1 =& G_{X_1}^{-1}(\zeta_1) \\
x_2 =& G_{X_2}^{-1}(\zeta_2)
\end{align}

### Exercise 2.3

Consider the input quantities $X_1,X_2$ with
* $X_1$ has best estimate 0.2 with uncertainty of 50%
* $X_2$ has best estimate -1.0 with uncertainty of 20%
* correlation between $X_1$ and $X_2$ is $\rho=-0.7$

Generate 1000 random draws using a bivariate normal copula function.


```python
# solution for exercise 2.3
from scipy.stats import norm

rho = -0.7    # correlation coefficient
draws = 1000  # number of random draws

# step 1: draw from bivariate normal distribution
z = multivariate_normal.rvs(zeros(2), array([[1.0, rho],[rho, 1.0]]), size=draws)

# step 2: evaluate cumulative distribution function of the copula
xi1 = norm.cdf(z[:,0])
xi2 = norm.cdf(z[:,1])

# step 3: evaluate inverse cumulative distr. function of the marginal distributions
x1 = norm.ppf(xi1, loc=0.2, scale = 0.1)
x2 = norm.ppf(xi2, loc=-1.0,scale = 0.2)

figure()
scatter(Xb[:,0],Xb[:,1],100,color="b", edgecolor="none",label="built-in function")
scatter(X[:,0],X[:,1],100, color="r", edgecolor="none",label="formula",alpha=0.2)
scatter(x1, x2, 100, color="g", edgecolor="none",label="copula",alpha=0.2)
legend();
```

### Exercise 2.4

Consider the input quantities $X_1, X_2$ with

* $X_1$ has best estimate $x_1=2.4$ with expanded uncertainty $U=0.4 (k=2)$ under normality assumption

* $X_2$ is located in $[-1.5, 1.5]$

* $X_1, X_2$ are correlated with $\rho = 0.4$

Draw 1000 samples from their joint probability distribution using a normal distribution copula function.


```python
# solution for exercise 2.4
from scipy.stats import norm, uniform, multivariate_normal
rho=0.4       # correlation coefficient
draws=1000    # number of random draws

# step 1: draw from bivariate normal distribution
z = multivariate_normal.rvs(zeros(2), array([[1.0, rho],[rho, 1.0]]), size=draws)

# step 2: evaluate cumulative distribution function of the copula
xi1 = norm.cdf(z[:,0])
xi2 = norm.cdf(z[:,1])

# step 3: evaluate inverse cumulative distr. function of the marginal distributions
x1 = norm.ppf(xi1, loc=2.4, scale=0.2)
x2 = uniform.ppf(xi2, loc=-1.5, scale=3)


figure()
subplot(221)
scatter(x1, x2, 100, edgecolor="none")
ylim(-2,2); xlim(1.6,3.2)
subplot(222)
hist(x2,bins=100, orientation="horizontal",edgecolor="none")
ylim(-2,2)
subplot(223)
hist(x1,bins=100,edgecolor="none")
xlim(1.6,3.2);
```

dc13ee1b648760b3f3d48ea5634a52c2eee28f97
ipynb
Jupyter Notebook
02 Drawing from multivariate distributions.ipynb
TanjaDorst/MonteCarloWorkshop
679a0387fea0a3f947af4e7cc0fafd59bee10d7b
MIT
02 Drawing from multivariate distributions.ipynb
TanjaDorst/MonteCarloWorkshop
679a0387fea0a3f947af4e7cc0fafd59bee10d7b
MIT
02 Drawing from multivariate distributions.ipynb
TanjaDorst/MonteCarloWorkshop
679a0387fea0a3f947af4e7cc0fafd59bee10d7b
MIT
Qwen/Qwen-72B
1. YES
2. YES