code: def rmv_normal_cov(mu, C, size=1):
      
      mu_size = np.shape(mu)
      if size == 1:
          return np.random.multivariate_normal(mu, C, size).reshape(mu_size)
      else:
          return np.random.multivariate_normal(
              mu, C, size).reshape((size,) + mu_size)
docs: Random multivariate normal variates.