Hello!

I use scipy.stats.sem a lot and I would love if it was able to take multiple axis arguments as many numpy functions can. Looking at the source code sem and other functions in scipy.stats are already implemented mostly in terms of numpy functions so it seems like it would only require changing some parts of the logic. Taking stats.sem as an example I think it would require changing
n = a.shape[axis]

to something like

n = product(a.shape[axis]) For masked arrays, a.count(axis) is used which already works with multiple axes.

Before I start work on a PR I wanted to ask if there is some reason that this change would be considered a bad idea. Otherwise, if I write a decent PR with tests, benchmarks and documentation that updates all the functions in scipy.stats that might reasonably take multiple axis arguments, is it likely to be accepted? Thanks, Tom