Replace `np.rollaxis()` with `np.moveaxis()` but it is not drop-in replacement

NumPy documentation [recommends](https://numpy.org/doc/stable/reference/generated/numpy.rollaxis.html#numpy.rollaxis) using [`np.moveaxis()`](https://numpy.org/doc/stable/reference/generated/numpy.moveaxis.html) instead of [`np.rollaxis()`](https://numpy.org/doc/stable/reference/generated/numpy.rollaxis.html). Currently `astropy` uses `np.rollaxis()` more than `np.moveaxis()`, but there is no reason we should not follow NumPy's recommendation. The lines of code that call `np.rollaxis()` can be found with `git grep 'rollaxis('`. The following sub-packages need to be updated:

- [ ] `modeling` – some work completed, but [needs review](https://github.com/astropy/astropy/pull/15727#issuecomment-1855094364)
- [ ] `stats` – #15727
- [ ] `time` – #15727

Pull requests can replace `np.rollaxis()` with `np.moveaxis()` in one or more sub-packages.

Start with `git grep 'rollaxis('` and inspect the remaining calls in the `modeling`, `stats`, and `time` sub-packages, along with the work and review discussion in PR #15727. Check each use against NumPy's `moveaxis` documentation and the relevant package tests; done means the targeted calls are updated without regressions.

Written by the indexing model from the issue text.