Random multivariate normal variates.

def rmv_normal_cov(mu, C, size=1):
    """
    Random multivariate normal variates.
    """
    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)