Shape issues with MvNormal.random method

[Sayam753]

This issue popped out in #4180 while writing tests for MixtureSameFamily distribution with MvNormal distribution. MvNormal.random method does not broadcast when sample_shape is equal to batch_shape and produces inconsistent shapes with different parametrizations.

Please provide a minimal, self-contained, and reproducible example.

>>> import numpy as np
>>> import pymc3 as pm
WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.
>>> mu = np.random.randn(10, 3)
>>> mat = np.random.randn(3, 3)
>>> cov = mat @ mat.T
>>> tau = np.linalg.inv(cov)
>>> chol = np.linalg.cholesky(cov)
>>> 
>>> comp_dists = pm.MvNormal.dist(mu=mu, chol=chol, shape=(10, 3))
>>> comp_dists.random(size=(10)).shape
(10, 3)
>>> comp_dists.random(size=(20)).shape
(20, 10, 3)
>>> comp_dists = pm.MvNormal.dist(mu=mu, cov=cov, shape=(10, 3))
>>> comp_dists.random(size=(10)).shape
(10, 3)
>>> comp_dists.random(size=(20)).shape
(20, 3)
>>> comp_dists = pm.MvNormal.dist(mu=mu, tau=tau, shape=(10, 3))
>>> comp_dists.random(size=(10)).shape
(10, 3)
>>> comp_dists.random(size=(20)).shape
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/sayam/.local/lib/python3.8/site-packages/pymc3/distributions/multivariate.py", line 303, in random
    return mu + transformed.T
ValueError: operands could not be broadcast together with shapes (10,3) (20,3) 

Please provide any additional information below.
On exploring, it turns out that there are known related issues #4075, #3706, #2848, #3829, #3758 . I think there needs to be major refactoring of MvNormal.random method.

Versions and main components

• PyMC3 Version: 3.9.3 (GitHub master)
• Theano Version: 1.0.9
• Python Version: 3.8.6
• Operating system: Ubuntu 18.04 LTS
• How did you install PyMC3: (conda/pip) pip

[AlexAndorra]

Thanks for reporting @Sayam753. There are indeed various issues with the MvNormal.random method and a major refactoring would certainly be very beneficial!