def METHOD_NAME(n: int) -> None:
    key_mean, key_sqrt = jr.split(_key, 2)
    mean = jr.uniform(key_mean, shape=(n,))
    sqrt = jr.uniform(key_sqrt, shape=(n, n))
    covariance = sqrt @ sqrt.T
    # check that cholesky does not error
    _L = jnp.linalg.cholesky(covariance)  # noqa: F841
    dist = GaussianDistribution(loc=mean, scale=cola.PSD(Dense(covariance)))
    assert approx_equal(dist.mean(), mean)
    assert approx_equal(dist.variance(), covariance.diagonal())
    assert approx_equal(dist.stddev(), jnp.sqrt(covariance.diagonal()))
    assert approx_equal(dist.covariance(), covariance)
    assert isinstance(dist.scale, Dense)
    assert cola.PSD in dist.scale.annotations
    y = jr.uniform(_key, shape=(n,))
    tfp_dist = MultivariateNormalFullCovariance(loc=mean, covariance_matrix=covariance)
    assert approx_equal(dist.log_prob(y), tfp_dist.log_prob(y))
    assert approx_equal(dist.kl_divergence(dist), 0.0)