[BUG] fix rvs output shape in multivariate distributions

[JDkuba]

Reference issue

Closes #7689

What does this implement/fix?

Some multivariate distributions (multivariate_normal, wishart and invwishart) rvs methods return output with incorrect shape. As @WarrenWeckesser noted, it is due to the indiscriminate use of squeeze.

Additional information

In my solution, axes that were given with the size parameter are being skipped by squeeze.

[rgommers]

CI failures are real:

>  out = out.squeeze(axis)
E       TypeError: function takes exactly 0 arguments (1 given)

[rgommers]

len(out.shape) --> out.ndim

[mdhaber]

We do want to fix this, but we need to be careful to warn users in advance. This PR would make the change immediately and it has merge conflicts, so I'll go ahead and close it. The need to fix this is tracked in gh-16278.