Bivariate normal distribution

I want to generate a bivariate normal distribution

\begin{equation}
  \begin{pmatrix}
    x \\
    y
  \end{pmatrix}
  = N
  \begin{pmatrix}
    \begin{pmatrix}
      \mu_x\\
      \mu_y
    \end{pmatrix}
    _,
    \begin{pmatrix}
      \sigma_x^2 & \rho\sigma_x\sigma_y \\
      \rho\sigma_x\sigma_y & \sigma_y^2 \\
    \end{pmatrix}
  \end{pmatrix}
\end{equation}

First, generate $ x = N (\ mu_x, \ sigma_x) $. Next, x should be fixed to generate y. The conditional distribution of y at this time is

y|x \sim N(\mu_y + \rho\frac{\sigma_y}{\sigma_x}(x-\mu_x), (1-\rho^2)\sigma_y^2))

Given in.

When implemented in python

python3

import numpy as np
import matplotlib.pyplot as plt

def MVNORM(mu_x, sigma_x, mu_y, sigma_y, rho, N=1):
  x = np.random.normal(mu_x, sigma_x, N)
  y = np.random.normal(mu_y + rho*sigma_y/sigma_x*(x-mu_x), np.sqrt((1-rho**2)*sigma_y**2), N)
  return([x,y])

#Mean 0, standard deviation 9 x and correlation coefficient 0.Generate 1000 y with standard deviation 3 with mean 3 correlated at 9
MV = MVNORM(0, 9, 3, 3, 0.9 ,1000)

plt.scatter(MV[0], MV[1])
plt.show()

SUMMARY = (np.mean(MV[0]), np.std(MV[0]), np.mean(MV[1]), np.std(MV[1]))
print("mean X: {0[0]:0.2f}, stdev X: {0[2]:0.2f}, mean Y: {0[1]:0.2f}, stdev Y: {0[3]:0.2f}".format(SUMMARY))

mean X: 0.14, stdev X: 9.09, mean Y: 2.99, stdev Y: 3.08