Improve MatrixNormal RandomVariable

[ricardoV94]

Description

This is the only random variable with ndim=2... and it's in need of some maintenance. It doesn't seem to adhere to the RandomVariable API regarding batched dimensions: #6556 (comment)

If mu has to be a matrix, it should be broadcasted exactly like rowchol and colchol, but it seems to follow it's own logic

pymc/pymc/distributions/multivariate.py

Lines 1653 to 1657 in c7279b5

	# Broadcasting all parameters
	shapes = [mu.shape, output_shape]
	broadcastable_shape = broadcast_dist_samples_shape(shapes, size=size)
	mu = np.broadcast_to(mu, shape=broadcastable_shape)
	rowchol = np.broadcast_to(rowchol, shape=size + rowchol.shape[-2:])

It also seems to not have enough tests regarding broadcasting behavior in the random number generation: #6556 (comment)

We should clean it up a bit. Perhaps it makes sense to move it temporarily to pymc-experimental?

Related to #6556