2021-11-29
2021-12-10

### Steps:

A widely used method for drawing (sampling) a random vector x from the N-dimensional multivariate normal distribution with mean vector μ and covariance matrix Σ works as follows:[35]

1. Find any real matrix A such that AAT = Σ. When Σ is positive-definite, the Cholesky decomposition is typically used, and the extended form of this decomposition can always be used (as the covariance matrix may be only positive semi-definite) in both cases a suitable matrix A is obtained. An alternative is to use the matrix A = ½ obtained from a spectral decomposition Σ = UΛU−1 of Σ. The former approach is more computationally straightforward but the matrices A change for different orderings of the elements of the random vector, while the latter approach gives matrices that are related by simple re-orderings. In theory both approaches give equally good ways of determining a suitable matrix A, but there are differences in computation time.
2. Let z = (z1, …, zN)T be a vector whose components are N independent standard normal variates (which can be generated, for example, by using the Box–Muller transform).
3. Let x be μ + Az. This has the desired distribution due to the affine transformation property.

Python Implementation:

1) From scratch

# # Define the desired distribution to sample from:

d = 2 # Number of dimensions
mean = np.matrix([[0.], [1.]])
covariance = np.matrix([
[1, 0.8],
[0.8, 1]
])

# # Compute the Decomposition:

A = np.linalg.cholesky(covariance)

# # Sample X from standard normal

n = 50 # Samples to draw
Z = np.random.normal(size=(d, n))

# # Apply the transformation

X = A.dot(Z) + mean

```# Plot the samples and the distribution
fig, ax = plt.subplots(figsize=(6, 4.5))
# Plot bivariate distribution
x1, x2, p = generate_surface(mean, covariance, d)
con = ax.contourf(x1, x2, p, 33, cmap=cm.YlGnBu)
# Plot samples
ax.plot(Y[0,:], Y[1,:], 'ro', alpha=.6,
markeredgecolor='k', markeredgewidth=0.5)
ax.set_xlabel('\$y_1\$', fontsize=13)
ax.set_ylabel('\$y_2\$', fontsize=13)
ax.axis([-2.5, 2.5, -1.5, 3.5])
ax.set_aspect('equal')
ax.set_title('Samples from bivariate normal distribution')
cbar = plt.colorbar(con)
cbar.ax.set_ylabel('density: \$p(y_1, y_2)\$', fontsize=13)
plt.show()

```

## 2) Using Numpy Sampler

Numpy has a built-in multivariate normal sampling function:

```z = np.random.multivariate_normal(mean=mean, cov=covariance, size=n)
```
```y = np.transpose(z)
```
```# Plot density function.
sns.jointplot(x=y[0], y=y[1], kind="kde", space=0);
```
##### Amir Masoud Sefidian
Data Scientist, Machine Learning Engineer, Researcher, Software Developer