itudes) transforms a set of
uncorrelated random variables into a new set of random variables with
the desired covariance. When a coloring transform is applied to a
sample of points distributed according to a multivariate normal
distribution with identity covariance and zero mean, the covariance of
the transformed sample is approximately the covariance matrix used
in the coloring transform.

Parameters
----------
x : array_like
    An array of points. The last dimension must correspond with the
    dimensionality of the space, i.e., the number of columns in the
    covariance matrix.

Returns
-------
x_ : array_like
    The transformed array of points.

References
----------
.. [1] "Whitening Transformation". Wikipedia.
       https://en.wikipedia.org/wiki/Whitening_transformation
.. [2] Novak, Lukas, and Miroslav Vorechovsky. "Generalization of
       coloring linear transformation". Transactions of VSB 18.2
       (2018): 31-35. :doi:`10.31490/tces-2018-0013`

Examples
--------
>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng(1638083107694713882823079058616272161)
>>> n = 3
>>> A = rng.random(size=(n, n))
>>> cov_array = A @ A.T  # make matrix symmetric positive definite
>>> cholesky = np.linalg.cholesky(cov_array)
>>> cov_object = stats.Covariance.from_cholesky(cholesky)
>>> x = rng.multivariate_normal(np.zeros(n), np.eye(n), size=(10000))
>>> x_ = cov_object.colorize(x)
>>> cov_data = np.cov(x_, rowvar=False)
>>> np.allclose(cov_data, cov_array, rtol=3e-2)
True
)