#### MIT License (c) 2018 by Andrew Lyasoff 

#### Jupyter notebook written in Python 3. It illustrates the simulation of multivariate normal samples with a given covariance matrix, and also the construction of univariate histograms.


```python
from sympy import init_printing
init_printing()
import json
import numpy as np
from numpy import *
```

Read the file that contains the covariance matrix (adjust the path in the next line accordingly).


```python
file = open('Cov.json','r')
Cov=array(json.load(file))
Cov
```




    array([[ 0.034969 ,  0.0314721,  0.       ,  0.       ],
           [ 0.0314721,  0.034969 ,  0.       ,  0.       ],
           [ 0.       ,  0.       ,  1.0630576, -0.5315288],
           [ 0.       ,  0.       , -0.5315288,  0.559504 ]])



First check if this is indeed a covariance matrix:


```python
Cov-Cov.T
```




    array([[0., 0., 0., 0.],
           [0., 0., 0., 0.],
           [0., 0., 0., 0.],
           [0., 0., 0., 0.]])



It is indeed selfadjoint. Now check its rank.


```python
rnk=np.linalg.matrix_rank(Cov)
rnk
```




    4



You can also find the eigenvalues and the associated eigenvectors.


```python
np.linalg.eig(Cov)
```




    (array([0.0664411, 0.0034969, 1.3994257, 0.2231359]),
     array([[ 0.70710678, -0.70710678,  0.        ,  0.        ],
            [ 0.70710678,  0.70710678,  0.        ,  0.        ],
            [ 0.        ,  0.        ,  0.84501075,  0.53474931],
            [ 0.        ,  0.        , -0.53474931,  0.84501075]]))



Or simply compute the determinant.


```python
np.linalg.det(Cov)
```

Now find the Cholesky decomposition:


```python
CovU=np.linalg.cholesky(Cov)
CovU
```




    array([[ 0.187     ,  0.        ,  0.        ,  0.        ],
           [ 0.1683    ,  0.08151141,  0.        ,  0.        ],
           [ 0.        ,  0.        ,  1.03104685,  0.        ],
           [ 0.        ,  0.        , -0.51552342,  0.54197749]])



Test if this is what you want:


```python
dot((CovU),CovU.T)-Cov
```




    array([[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.11022302e-16]])



The next function generates a Gaussian vector with independent standard normal (${\cal N}(0,1)$) components (the exact dimension is supplied as an argument).


```python
def gauss_vec(dimen):
    vec=[]
    for i in range(dimen):
        vec+=[np.random.normal(0,1)]
    return vec

```


```python
gauss_vec(3)
```


```python
gauss_vec(rnk)
```

Now create a sample from a Gaussian vector with covariance matrix Cov.


```python
dot((CovU),gauss_vec(rnk))
```




    array([ 0.09825125,  0.06657924, -0.99930286,  0.47786305])




```python
list(dot((CovU),gauss_vec(rnk)))
```

Now we do the same by using the $LDL^\dagger$ decomposition.


```python
D,U=np.linalg.eig(Cov)
```


```python
CovL=dot(U,diag(sqrt(D)))
```

Check if this indeed the intended matrix:


```python
dot((CovL),CovL.T)-Cov
```




    array([[-6.93889390e-18,  6.93889390e-18,  0.00000000e+00,
             0.00000000e+00],
           [ 6.93889390e-18,  6.93889390e-18,  0.00000000e+00,
             0.00000000e+00],
           [ 0.00000000e+00,  0.00000000e+00,  2.22044605e-16,
            -1.11022302e-16],
           [ 0.00000000e+00,  0.00000000e+00, -1.11022302e-16,
             1.11022302e-16]])




```python
list(dot((Cov),gauss_vec(rnk)))
```


```python

```

531b6eb58027680246ca34b2608d7a1285fa171f
ipynb
Jupyter Notebook
Multivariate_Normal_Dist_Examples_Python.ipynb
AndrewLyasoff/SMAP
6eeea8953a26a05b1e23387067109d23b2011824
MIT
2018-09-04T19:12:32.000Z
2022-03-20T02:05:44.000Z
Multivariate_Normal_Dist_Examples_Python.ipynb
lhyzh/SMAP
f6687291769d4c16a0d51a06a941384f646bb432
MIT
2019-02-13T14:24:38.000Z
2019-02-13T14:24:38.000Z
Multivariate_Normal_Dist_Examples_Python.ipynb
lhyzh/SMAP
f6687291769d4c16a0d51a06a941384f646bb432
MIT
2019-02-10T03:43:40.000Z
2021-03-28T03:53:22.000Z
Qwen/Qwen-72B
1. YES
2. YES