diff --git a/scipy-stubs/linalg/_matfuncs.pyi b/scipy-stubs/linalg/_matfuncs.pyi
index 73ce35c5..c9b356a8 100644
--- a/scipy-stubs/linalg/_matfuncs.pyi
+++ b/scipy-stubs/linalg/_matfuncs.pyi
@@ -3,7 +3,6 @@ from typing import Any, Literal, TypeAlias, overload
 
 import numpy as np
 import optype.numpy as onp
-import optype.typing as opt
 from ._expm_frechet import expm_cond, expm_frechet
 from ._matfuncs_sqrtm import sqrtm
 
@@ -52,7 +51,7 @@ def fractional_matrix_power(A: onp.ToFloat2D, t: onp.ToInt) -> _Real2D: ...
 @overload  # complex, int
 def fractional_matrix_power(A: onp.ToComplex2D, t: onp.ToInt) -> _Numeric2D: ...
 @overload  # complex, float
-def fractional_matrix_power(A: onp.ToComplex2D, t: opt.Just[float] | np.floating[Any]) -> _Complex2D: ...
+def fractional_matrix_power(A: onp.ToComplex2D, t: onp.ToJustFloat) -> _Complex2D: ...
 
 # NOTE: return dtype depends on the sign of the values
 @overload  # disp: True = ...
diff --git a/scipy-stubs/linalg/_sketches.pyi b/scipy-stubs/linalg/_sketches.pyi
index 28279bf4..561c585c 100644
--- a/scipy-stubs/linalg/_sketches.pyi
+++ b/scipy-stubs/linalg/_sketches.pyi
@@ -1,39 +1,28 @@
-from collections.abc import Sequence
-from typing import Any, TypeAlias, TypeVar, overload
+from typing import overload
 
 import numpy as np
 import optype.numpy as onp
-import optype.typing as opt
-import scipy._typing as spt
+from scipy._typing import ToRNG
 from scipy.sparse import csc_matrix
 
 __all__ = ["clarkson_woodruff_transform"]
 
-_ST = TypeVar("_ST", bound=np.generic)
-_VT = TypeVar("_VT")
-
-_ToJust2D: TypeAlias = onp.CanArrayND[_ST] | Sequence[onp.CanArrayND[_ST]] | Sequence[Sequence[opt.Just[_VT] | _ST]]
-
 ###
 
-def cwt_matrix(n_rows: onp.ToInt, n_columns: onp.ToInt, rng: spt.ToRNG = None) -> csc_matrix: ...
+def cwt_matrix(n_rows: onp.ToInt, n_columns: onp.ToInt, rng: ToRNG = None) -> csc_matrix: ...
 
 #
 @overload
-def clarkson_woodruff_transform(
-    input_matrix: onp.ToInt2D,
-    sketch_size: onp.ToInt,
-    rng: spt.ToRNG = None,
-) -> onp.Array2D[np.int_]: ...
+def clarkson_woodruff_transform(input_matrix: onp.ToInt2D, sketch_size: onp.ToInt, rng: ToRNG = None) -> onp.Array2D[np.int_]: ...
 @overload
 def clarkson_woodruff_transform(
-    input_matrix: _ToJust2D[np.floating[Any], float],
+    input_matrix: onp.ToJustFloat2D,
     sketch_size: onp.ToInt,
-    rng: spt.ToRNG = None,
+    rng: ToRNG = None,
 ) -> onp.Array2D[np.float64 | np.longdouble]: ...
 @overload
 def clarkson_woodruff_transform(
-    input_matrix: _ToJust2D[np.complexfloating[Any, Any], complex],
+    input_matrix: onp.ToJustComplex2D,
     sketch_size: onp.ToInt,
-    rng: spt.ToRNG = None,
+    rng: ToRNG = None,
 ) -> onp.Array2D[np.complex64 | np.clongdouble]: ...
diff --git a/scipy-stubs/linalg/_special_matrices.pyi b/scipy-stubs/linalg/_special_matrices.pyi
index 0b1071f8..77703031 100644
--- a/scipy-stubs/linalg/_special_matrices.pyi
+++ b/scipy-stubs/linalg/_special_matrices.pyi
@@ -1,10 +1,11 @@
 from collections.abc import Sequence
-from typing import Any, Literal as L, TypeAlias, overload
+from typing import Literal as L, TypeAlias, overload
 from typing_extensions import TypeVar, deprecated
 
 import numpy as np
 import numpy.typing as npt
 import optype.numpy as onp
+import optype.numpy.compat as npc
 import optype.typing as opt
 from scipy._typing import ConvMode
 
@@ -28,273 +29,270 @@ __all__ = [
     "toeplitz",
 ]
 
-_SCT = TypeVar("_SCT", bound=np.generic, default=np.number[Any] | np.bool_ | np.object_)
+_SCT = TypeVar("_SCT", bound=np.generic, default=npc.number | np.bool_ | np.object_)
 
 _Kind: TypeAlias = L["symmetric", "upper", "lower"]
 
-_Array2d: TypeAlias = onp.Array2D[_SCT]
-_Array2dPlus: TypeAlias = onp.Array[onp.AtLeast2D, _SCT]
-_Array3dPlus: TypeAlias = onp.Array[onp.AtLeast3D, _SCT]
+_Array2ND: TypeAlias = onp.Array[onp.AtLeast2D, _SCT]
+_Array3ND: TypeAlias = onp.Array[onp.AtLeast3D, _SCT]
 
-_Int2d: TypeAlias = _Array2d[np.int_]
-_Float2d: TypeAlias = _Array2d[np.float64]
-_Complex2d: TypeAlias = _Array2d[np.complex128]
+_Int2D: TypeAlias = onp.Array2D[np.int_]
+_Int3ND: TypeAlias = _Array3ND[np.int_]
+_Float2D: TypeAlias = onp.Array2D[np.float64]
+_Float3ND: TypeAlias = _Array3ND[np.float64]
+_Complex2D: TypeAlias = onp.Array2D[np.complex128]
+_Complex3ND: TypeAlias = _Array3ND[np.complex128]
 
-_Int3dPlus: TypeAlias = _Array3dPlus[np.int_]
-_Float3dPlus: TypeAlias = _Array3dPlus[np.float64]
-_Complex3dPlus: TypeAlias = _Array3dPlus[np.complex128]
+_To0D: TypeAlias = _SCT | onp.CanArray0D[_SCT]
+_To1D: TypeAlias = Sequence[_To0D[_SCT]] | onp.CanArrayND[_SCT]
+_ToStrict1D: TypeAlias = Sequence[_To0D[_SCT]] | onp.CanArray1D[_SCT]
+_ToStrict2ND: TypeAlias = onp.SequenceND[_To1D[_SCT]] | onp.CanArrayND[_SCT, onp.AtLeast2D]
+_ToND: TypeAlias = onp.SequenceND[_To0D[_SCT]] | onp.SequenceND[_To1D[_SCT]] | onp.CanArrayND[_SCT]
+_ToDType: TypeAlias = type[_SCT] | np.dtype[_SCT] | onp.HasDType[np.dtype[_SCT]]
 
-_ToArray1d: TypeAlias = onp.CanArrayND[_SCT] | Sequence[_SCT]
-_ToArrayNd: TypeAlias = onp.CanArrayND[_SCT] | onp.SequenceND[onp.CanArrayND[_SCT]] | onp.SequenceND[_SCT]
-_ToArrayStrict1d: TypeAlias = onp.CanArray1D[_SCT] | Sequence[_SCT]
-_ToArrayStrict1dPlus: TypeAlias = onp.CanArrayND[_SCT, onp.AtLeast2D] | onp.SequenceND[Sequence[_SCT] | onp.CanArrayND[_SCT]]
-
-_JustInt1d: TypeAlias = Sequence[opt.JustInt | np.int_]
-_JustFloat1d: TypeAlias = Sequence[opt.Just[float] | np.float64]
-_JustComplex1d: TypeAlias = Sequence[opt.Just[complex] | np.complex128]
-_JustInt2dPlus: TypeAlias = onp.SequenceND[_JustInt1d]
-_JustFloat2dPlus: TypeAlias = onp.SequenceND[_JustFloat1d]
-_JustComplex2dPlus: TypeAlias = onp.SequenceND[_JustComplex1d]
+_ToJustIntStrict2ND: TypeAlias = onp.SequenceND[onp.ToJustInt1D] | onp.CanArrayND[npc.integer, onp.AtLeast2D]
+_ToJustFloatStrict2ND: TypeAlias = onp.SequenceND[onp.ToJustFloat1D] | onp.CanArrayND[npc.floating, onp.AtLeast2D]
+_ToJustComplexStrict2ND: TypeAlias = onp.SequenceND[onp.ToJustComplex1D] | onp.CanArrayND[npc.complexfloating, onp.AtLeast2D]
 
 ###
 
 #
 @overload
 @deprecated("`kron` has been deprecated in favour of `numpy.kron` in SciPy 1.15.0 and will be removed in SciPy 1.17.0.")
-def kron(a: _Array2d[_SCT], b: _Array2d[_SCT]) -> _Array2d[_SCT]: ...
+def kron(a: onp.Array2D[_SCT], b: onp.Array2D[_SCT]) -> onp.Array2D[_SCT]: ...
 @overload
 @deprecated("`kron` has been deprecated in favour of `numpy.kron` in SciPy 1.15.0 and will be removed in SciPy 1.17.0.")
 def kron(a: onp.ArrayND[_SCT], b: onp.ArrayND[_SCT]) -> onp.Array[onp.AtLeast2D, _SCT]: ...
 
 #
 @overload
-def toeplitz(c: _JustInt1d, r: _JustInt1d | None = None) -> _Int2d: ...
+def toeplitz(c: onp.ToJustInt1D, r: onp.ToJustInt1D | None = None) -> _Int2D: ...
 @overload
-def toeplitz(c: _JustFloat1d, r: _JustFloat1d | None = None) -> _Float2d: ...
+def toeplitz(c: onp.ToJustFloat1D, r: onp.ToJustFloat1D | None = None) -> _Float2D: ...
 @overload
-def toeplitz(c: _JustComplex1d, r: _JustComplex1d | None = None) -> _Complex2d: ...
+def toeplitz(c: onp.ToJustComplex1D, r: onp.ToJustComplex1D | None = None) -> _Complex2D: ...
 @overload
-def toeplitz(c: _ToArrayStrict1d[_SCT], r: _ToArrayStrict1d[_SCT] | None = None) -> _Array2d[_SCT]: ...
+def toeplitz(c: _ToStrict1D[_SCT], r: _ToStrict1D[_SCT] | None = None) -> onp.Array2D[_SCT]: ...
 @overload
 @deprecated(
     "Beginning in SciPy 1.17, multidimensional input will be treated as a batch, not `ravel`ed. "
     "To preserve the existing behavior and silence this warning, `ravel` aruments before passing them to `toeplitz`.",
     category=FutureWarning,
 )
-def toeplitz(c: _JustInt2dPlus, r: _JustInt2dPlus | None = None) -> _Int2d: ...
+def toeplitz(c: _ToJustIntStrict2ND, r: _ToJustIntStrict2ND | None = None) -> _Int2D: ...
 @overload
 @deprecated(
     "Beginning in SciPy 1.17, multidimensional input will be treated as a batch, not `ravel`ed. "
     "To preserve the existing behavior and silence this warning, `ravel` aruments before passing them to `toeplitz`.",
     category=FutureWarning,
 )
-def toeplitz(c: _JustFloat2dPlus, r: _JustFloat2dPlus | None = None) -> _Float2d: ...
+def toeplitz(c: _ToJustFloatStrict2ND, r: _ToJustFloatStrict2ND | None = None) -> _Float2D: ...
 @overload
 @deprecated(
     "Beginning in SciPy 1.17, multidimensional input will be treated as a batch, not `ravel`ed. "
     "To preserve the existing behavior and silence this warning, `ravel` aruments before passing them to `toeplitz`.",
     category=FutureWarning,
 )
-def toeplitz(c: _JustComplex2dPlus, r: _JustComplex2dPlus | None = None) -> _Complex2d: ...
+def toeplitz(c: _ToJustComplexStrict2ND, r: _ToJustComplexStrict2ND | None = None) -> _Complex2D: ...
 @overload
 @deprecated(
     "Beginning in SciPy 1.17, multidimensional input will be treated as a batch, not `ravel`ed. "
     "To preserve the existing behavior and silence this warning, `ravel` aruments before passing them to `toeplitz`.",
     category=FutureWarning,
 )
-def toeplitz(c: _ToArrayStrict1dPlus[_SCT], r: _ToArrayStrict1dPlus[_SCT] | None = None) -> _Array2d[_SCT]: ...
+def toeplitz(c: _ToStrict2ND[_SCT], r: _ToStrict2ND[_SCT] | None = None) -> onp.Array2D[_SCT]: ...
 @overload
-def toeplitz(c: _ToArray1d[_SCT], r: _ToArray1d[_SCT] | None = None) -> _Array2d[_SCT]: ...
+def toeplitz(c: _To1D[_SCT], r: _To1D[_SCT] | None = None) -> onp.Array2D[_SCT]: ...
 
 #
 @overload
-def circulant(c: _JustInt1d) -> _Int2d: ...
+def circulant(c: onp.ToJustIntStrict1D) -> _Int2D: ...
 @overload
-def circulant(c: _JustInt2dPlus) -> _Int3dPlus: ...
+def circulant(c: _ToJustIntStrict2ND) -> _Int3ND: ...
 @overload
-def circulant(c: _JustFloat1d) -> _Float2d: ...
+def circulant(c: onp.ToJustFloatStrict1D) -> _Float2D: ...
 @overload
-def circulant(c: _JustFloat2dPlus) -> _Float3dPlus: ...
+def circulant(c: _ToJustFloatStrict2ND) -> _Float3ND: ...
 @overload
-def circulant(c: _JustComplex1d) -> _Complex2d: ...
+def circulant(c: onp.ToJustComplexStrict1D) -> _Complex2D: ...
 @overload
-def circulant(c: _JustComplex2dPlus) -> _Complex3dPlus: ...
+def circulant(c: _ToJustComplexStrict2ND) -> _Complex3ND: ...
 @overload
-def circulant(c: _ToArrayStrict1d[_SCT]) -> _Array2d[_SCT]: ...
+def circulant(c: _ToStrict1D[_SCT]) -> onp.Array2D[_SCT]: ...
 @overload
-def circulant(c: _ToArrayStrict1dPlus[_SCT]) -> _Array3dPlus[_SCT]: ...
+def circulant(c: _ToStrict2ND[_SCT]) -> _Array3ND[_SCT]: ...
 @overload
-def circulant(c: _ToArrayNd[_SCT]) -> _Array2dPlus[_SCT]: ...
+def circulant(c: _ToND[_SCT]) -> _Array2ND[_SCT]: ...
 
 #
 @overload
-def companion(a: _JustInt1d) -> _Int2d: ...
+def companion(a: onp.ToJustIntStrict1D) -> _Int2D: ...
 @overload
-def companion(a: _JustInt2dPlus) -> _Int3dPlus: ...
+def companion(a: _ToJustIntStrict2ND) -> _Int3ND: ...
 @overload
-def companion(a: _JustFloat1d) -> _Float2d: ...
+def companion(a: onp.ToJustFloatStrict1D) -> _Float2D: ...
 @overload
-def companion(a: _JustFloat2dPlus) -> _Float3dPlus: ...
+def companion(a: _ToJustFloatStrict2ND) -> _Float3ND: ...
 @overload
-def companion(a: _JustComplex1d) -> _Complex2d: ...
+def companion(a: onp.ToJustComplexStrict1D) -> _Complex2D: ...
 @overload
-def companion(a: _JustComplex2dPlus) -> _Complex3dPlus: ...
+def companion(a: _ToJustComplexStrict2ND) -> _Complex3ND: ...
 @overload
-def companion(a: _ToArrayStrict1d[_SCT]) -> _Array2d[_SCT]: ...
+def companion(a: _ToStrict1D[_SCT]) -> onp.Array2D[_SCT]: ...
 @overload
-def companion(a: _ToArrayStrict1dPlus[_SCT]) -> _Array3dPlus[_SCT]: ...
+def companion(a: _ToStrict2ND[_SCT]) -> _Array3ND[_SCT]: ...
 @overload
-def companion(a: _ToArrayNd[_SCT]) -> _Array2dPlus[_SCT]: ...
+def companion(a: _ToND[_SCT]) -> _Array2ND[_SCT]: ...
 
 #
 @overload
-def convolution_matrix(a: _JustInt1d, n: onp.ToInt, mode: ConvMode = "full") -> _Int2d: ...
+def convolution_matrix(a: onp.ToJustIntStrict1D, n: onp.ToInt, mode: ConvMode = "full") -> _Int2D: ...
 @overload
-def convolution_matrix(a: _JustInt2dPlus, n: onp.ToInt, mode: ConvMode = "full") -> _Int3dPlus: ...
+def convolution_matrix(a: _ToJustIntStrict2ND, n: onp.ToInt, mode: ConvMode = "full") -> _Int3ND: ...
 @overload
-def convolution_matrix(a: _JustFloat1d, n: onp.ToInt, mode: ConvMode = "full") -> _Float2d: ...
+def convolution_matrix(a: onp.ToJustFloatStrict1D, n: onp.ToInt, mode: ConvMode = "full") -> _Float2D: ...
 @overload
-def convolution_matrix(a: _JustFloat2dPlus, n: onp.ToInt, mode: ConvMode = "full") -> _Float3dPlus: ...
+def convolution_matrix(a: _ToJustFloatStrict2ND, n: onp.ToInt, mode: ConvMode = "full") -> _Float3ND: ...
 @overload
-def convolution_matrix(a: _JustComplex1d, n: onp.ToInt, mode: ConvMode = "full") -> _Complex2d: ...
+def convolution_matrix(a: onp.ToJustComplexStrict1D, n: onp.ToInt, mode: ConvMode = "full") -> _Complex2D: ...
 @overload
-def convolution_matrix(a: _JustComplex2dPlus, n: onp.ToInt, mode: ConvMode = "full") -> _Complex3dPlus: ...
+def convolution_matrix(a: _ToJustComplexStrict2ND, n: onp.ToInt, mode: ConvMode = "full") -> _Complex3ND: ...
 @overload
-def convolution_matrix(a: _ToArrayStrict1d[_SCT], n: onp.ToInt, mode: ConvMode = "full") -> _Array2d[_SCT]: ...
+def convolution_matrix(a: _ToStrict1D[_SCT], n: onp.ToInt, mode: ConvMode = "full") -> onp.Array2D[_SCT]: ...
 @overload
-def convolution_matrix(a: _ToArrayStrict1dPlus[_SCT], n: onp.ToInt, mode: ConvMode = "full") -> _Array3dPlus[_SCT]: ...
+def convolution_matrix(a: _ToStrict2ND[_SCT], n: onp.ToInt, mode: ConvMode = "full") -> _Array3ND[_SCT]: ...
 @overload
-def convolution_matrix(a: _ToArrayNd[_SCT], n: onp.ToInt, mode: ConvMode = "full") -> _Array2dPlus[_SCT]: ...
+def convolution_matrix(a: _ToND[_SCT], n: onp.ToInt, mode: ConvMode = "full") -> _Array2ND[_SCT]: ...
 
 #
 @overload
-def fiedler(a: _JustInt1d) -> _Int2d: ...
+def fiedler(a: onp.ToJustIntStrict1D) -> _Int2D: ...
 @overload
-def fiedler(a: _JustInt2dPlus) -> _Int3dPlus: ...
+def fiedler(a: _ToJustIntStrict2ND) -> _Int3ND: ...
 @overload
-def fiedler(a: _JustFloat1d) -> _Float2d: ...
+def fiedler(a: onp.ToJustFloatStrict1D) -> _Float2D: ...
 @overload
-def fiedler(a: _JustFloat2dPlus) -> _Float3dPlus: ...
+def fiedler(a: _ToJustFloatStrict2ND) -> _Float3ND: ...
 @overload
-def fiedler(a: _JustComplex1d) -> _Complex2d: ...
+def fiedler(a: onp.ToJustComplexStrict1D) -> _Complex2D: ...
 @overload
-def fiedler(a: _JustComplex2dPlus) -> _Complex3dPlus: ...
+def fiedler(a: _ToJustComplexStrict2ND) -> _Complex3ND: ...
 @overload
-def fiedler(a: _ToArrayStrict1d[_SCT]) -> _Array2d[_SCT]: ...
+def fiedler(a: _ToStrict1D[_SCT]) -> onp.Array2D[_SCT]: ...
 @overload
-def fiedler(a: _ToArrayStrict1dPlus[_SCT]) -> _Array3dPlus[_SCT]: ...
+def fiedler(a: _ToStrict2ND[_SCT]) -> _Array3ND[_SCT]: ...
 @overload
-def fiedler(a: _ToArrayNd[_SCT]) -> _Array2dPlus[_SCT]: ...
+def fiedler(a: _ToND[_SCT]) -> _Array2ND[_SCT]: ...
 
 #
 @overload
-def fiedler_companion(a: _JustInt1d) -> _Int2d: ...
+def fiedler_companion(a: onp.ToJustIntStrict1D) -> _Int2D: ...
 @overload
-def fiedler_companion(a: _JustInt2dPlus) -> _Int3dPlus: ...
+def fiedler_companion(a: _ToJustIntStrict2ND) -> _Int3ND: ...
 @overload
-def fiedler_companion(a: _JustFloat1d) -> _Float2d: ...
+def fiedler_companion(a: onp.ToJustFloatStrict1D) -> _Float2D: ...
 @overload
-def fiedler_companion(a: _JustFloat2dPlus) -> _Float3dPlus: ...
+def fiedler_companion(a: _ToJustFloatStrict2ND) -> _Float3ND: ...
 @overload
-def fiedler_companion(a: _JustComplex1d) -> _Complex2d: ...
+def fiedler_companion(a: onp.ToJustComplexStrict1D) -> _Complex2D: ...
 @overload
-def fiedler_companion(a: _JustComplex2dPlus) -> _Complex3dPlus: ...
+def fiedler_companion(a: _ToJustComplexStrict2ND) -> _Complex3ND: ...
 @overload
-def fiedler_companion(a: _ToArrayStrict1d[_SCT]) -> _Array2d[_SCT]: ...
+def fiedler_companion(a: _ToStrict1D[_SCT]) -> onp.Array2D[_SCT]: ...
 @overload
-def fiedler_companion(a: _ToArrayStrict1dPlus[_SCT]) -> _Array3dPlus[_SCT]: ...
+def fiedler_companion(a: _ToStrict2ND[_SCT]) -> _Array3ND[_SCT]: ...
 @overload
-def fiedler_companion(a: _ToArrayNd[_SCT]) -> _Array2dPlus[_SCT]: ...
+def fiedler_companion(a: _ToND[_SCT]) -> _Array2ND[_SCT]: ...
 
 # TODO(jorenham): batching
 @overload
-def leslie(f: _JustInt1d, s: _JustInt1d) -> _Int2d: ...
+def leslie(f: onp.ToJustIntStrict1D, s: onp.ToJustIntStrict1D) -> _Int2D: ...
 @overload
-def leslie(f: _JustInt1d, s: _JustInt2dPlus) -> _Int3dPlus: ...
+def leslie(f: onp.ToJustIntND, s: _ToJustIntStrict2ND) -> _Int3ND: ...
 @overload
-def leslie(f: _JustInt2dPlus, s: _JustInt1d | _JustInt2dPlus) -> _Int3dPlus: ...
+def leslie(f: _ToJustIntStrict2ND, s: onp.ToJustIntND) -> _Int3ND: ...
 @overload
-def leslie(f: _JustFloat1d, s: _JustFloat1d) -> _Float2d: ...
+def leslie(f: onp.ToJustFloatStrict1D, s: onp.ToJustFloatStrict1D) -> _Float2D: ...
 @overload
-def leslie(f: _JustFloat1d, s: _JustFloat2dPlus) -> _Float3dPlus: ...
+def leslie(f: onp.ToJustFloatND, s: _ToJustFloatStrict2ND) -> _Float3ND: ...
 @overload
-def leslie(f: _JustFloat2dPlus, s: _JustFloat1d | _JustFloat2dPlus) -> _Float3dPlus: ...
+def leslie(f: _ToJustFloatStrict2ND, s: onp.ToJustFloatND) -> _Float3ND: ...
 @overload
-def leslie(f: _JustComplex1d, s: _JustComplex1d) -> _Complex2d: ...
+def leslie(f: onp.ToJustComplexStrict1D, s: onp.ToJustComplexStrict1D) -> _Complex2D: ...
 @overload
-def leslie(f: _JustComplex1d, s: _JustComplex2dPlus) -> _Complex3dPlus: ...
+def leslie(f: onp.ToJustComplexND, s: _ToJustComplexStrict2ND) -> _Complex3ND: ...
 @overload
-def leslie(f: _JustComplex2dPlus, s: _JustComplex1d | _JustComplex2dPlus) -> _Complex3dPlus: ...
+def leslie(f: _ToJustComplexStrict2ND, s: onp.ToJustComplexND) -> _Complex3ND: ...
 @overload
-def leslie(f: _ToArrayStrict1d[_SCT], s: _ToArrayStrict1d[_SCT]) -> _Array2d[_SCT]: ...
+def leslie(f: _ToStrict1D[_SCT], s: _ToStrict1D[_SCT]) -> onp.Array2D[_SCT]: ...
 @overload
-def leslie(f: _ToArrayStrict1d[_SCT], s: _ToArrayStrict1dPlus[_SCT]) -> _Array3dPlus[_SCT]: ...
+def leslie(f: _ToStrict1D[_SCT], s: _ToStrict2ND[_SCT]) -> _Array3ND[_SCT]: ...
 @overload
-def leslie(f: _ToArrayStrict1dPlus[_SCT], s: _ToArrayNd[_SCT]) -> _Array3dPlus[_SCT]: ...
+def leslie(f: _ToStrict2ND[_SCT], s: _ToND[_SCT]) -> _Array3ND[_SCT]: ...
 @overload
-def leslie(f: _ToArrayNd[_SCT], s: _ToArrayNd[_SCT]) -> _Array2dPlus[_SCT]: ...
+def leslie(f: _ToND[_SCT], s: _ToND[_SCT]) -> _Array2ND[_SCT]: ...
 
 #
 @overload
-def block_diag() -> _Float2d: ...  # shape=(1, 0)
+def block_diag() -> _Float2D: ...  # shape=(1, 0)
 @overload
-def block_diag(arr0: _JustInt1d, /, *arrs: _JustInt1d) -> _Int2d: ...
+def block_diag(arr0: onp.ToJustInt1D, /, *arrs: onp.ToJustInt1D) -> _Int2D: ...
 @overload
-def block_diag(arr0: _JustFloat1d, /, *arrs: _JustFloat1d) -> _Float2d: ...
+def block_diag(arr0: onp.ToJustFloat1D, /, *arrs: onp.ToJustFloat1D) -> _Float2D: ...
 @overload
-def block_diag(arr0: _JustComplex1d, /, *arrs: _JustComplex1d) -> _Complex2d: ...
+def block_diag(arr0: onp.ToJustComplex1D, /, *arrs: onp.ToJustComplex1D) -> _Complex2D: ...
 @overload
-def block_diag(arr0: _ToArray1d[_SCT], /, *arrs: _ToArray1d[_SCT]) -> _Array2d[_SCT]: ...
+def block_diag(arr0: _To1D[_SCT], /, *arrs: _To1D[_SCT]) -> onp.Array2D[_SCT]: ...
 
 #
-def dft(n: onp.ToInt, scale: L["sqrtn", "n"] | None = None) -> _Complex2d: ...
+def dft(n: onp.ToInt, scale: L["sqrtn", "n"] | None = None) -> _Complex2D: ...
 
 #
 @overload
-def hadamard(n: onp.ToInt, dtype: type[opt.JustInt]) -> _Int2d: ...
+def hadamard(n: onp.ToInt, dtype: type[opt.JustInt]) -> _Int2D: ...
 @overload
-def hadamard(n: onp.ToInt, dtype: type[opt.Just[float]]) -> _Float2d: ...
+def hadamard(n: onp.ToInt, dtype: type[opt.JustFloat]) -> _Float2D: ...
 @overload
-def hadamard(n: onp.ToInt, dtype: type[opt.Just[complex]]) -> _Complex2d: ...
+def hadamard(n: onp.ToInt, dtype: type[opt.JustComplex]) -> _Complex2D: ...
 @overload
-def hadamard(n: onp.ToInt, dtype: type[_SCT] | np.dtype[_SCT] | onp.HasDType[np.dtype[_SCT]]) -> _Array2d[_SCT]: ...
+def hadamard(n: onp.ToInt, dtype: _ToDType[_SCT]) -> onp.Array2D[_SCT]: ...
 @overload
-def hadamard(n: onp.ToInt, dtype: npt.DTypeLike = ...) -> _Array2d[np.generic]: ...
+def hadamard(n: onp.ToInt, dtype: npt.DTypeLike = ...) -> onp.Array2D: ...
 
 #
 @overload
-def hankel(c: _JustInt1d, r: _JustInt1d | None = None) -> _Int2d: ...
+def hankel(c: onp.ToJustInt1D, r: onp.ToJustInt1D | None = None) -> _Int2D: ...
 @overload
-def hankel(c: _JustFloat1d, r: _JustFloat1d | None = None) -> _Float2d: ...
+def hankel(c: onp.ToJustFloat1D, r: onp.ToJustFloat1D | None = None) -> _Float2D: ...
 @overload
-def hankel(c: _JustComplex1d, r: _JustComplex1d | None = None) -> _Complex2d: ...
+def hankel(c: onp.ToJustComplex1D, r: onp.ToJustComplex1D | None = None) -> _Complex2D: ...
 @overload
-def hankel(c: _ToArray1d[_SCT], r: _ToArray1d[_SCT] | None = None) -> _Array2d[_SCT]: ...
+def hankel(c: _To1D[_SCT], r: _To1D[_SCT] | None = None) -> onp.Array2D[_SCT]: ...
 
 #
-def helmert(n: onp.ToInt, full: bool = False) -> _Float2d: ...
+def helmert(n: onp.ToInt, full: bool = False) -> _Float2D: ...
 
 #
-def hilbert(n: onp.ToInt) -> _Float2d: ...
+def hilbert(n: onp.ToInt) -> _Float2D: ...
 
 #
 @overload
-def invhilbert(n: onp.ToInt, exact: L[False] = False) -> _Float2d: ...
+def invhilbert(n: onp.ToInt, exact: L[False] = False) -> _Float2D: ...
 @overload
-def invhilbert(n: onp.ToInt, exact: L[True]) -> _Array2d[np.int64] | _Array2d[np.object_]: ...
+def invhilbert(n: onp.ToInt, exact: L[True]) -> onp.Array2D[np.int64 | np.object_]: ...
 
 #
 @overload
-def pascal(n: onp.ToInt, kind: _Kind = "symmetric", exact: L[True] = True) -> _Array2d[np.int64] | _Array2d[np.object_]: ...
+def pascal(n: onp.ToInt, kind: _Kind = "symmetric", exact: L[True] = True) -> onp.Array2D[np.int64 | np.object_]: ...
 @overload
-def pascal(n: onp.ToInt, kind: _Kind, exact: L[False]) -> _Float2d: ...
+def pascal(n: onp.ToInt, kind: _Kind, exact: L[False]) -> _Float2D: ...
 @overload
-def pascal(n: onp.ToInt, kind: _Kind = "symmetric", *, exact: L[False]) -> _Float2d: ...
+def pascal(n: onp.ToInt, kind: _Kind = "symmetric", *, exact: L[False]) -> _Float2D: ...
 
 #
 @overload
-def invpascal(n: onp.ToInt, kind: _Kind = "symmetric", exact: L[True] = True) -> _Array2d[np.int64] | _Array2d[np.object_]: ...
+def invpascal(n: onp.ToInt, kind: _Kind = "symmetric", exact: L[True] = True) -> onp.Array2D[np.int64 | np.object_]: ...
 @overload
-def invpascal(n: onp.ToInt, kind: _Kind, exact: L[False]) -> _Float2d: ...
+def invpascal(n: onp.ToInt, kind: _Kind, exact: L[False]) -> _Float2D: ...
 @overload
-def invpascal(n: onp.ToInt, kind: _Kind = "symmetric", *, exact: L[False]) -> _Float2d: ...
+def invpascal(n: onp.ToInt, kind: _Kind = "symmetric", *, exact: L[False]) -> _Float2D: ...
diff --git a/scipy-stubs/sparse/_base.pyi b/scipy-stubs/sparse/_base.pyi
index 0f7b9131..1e0d04cb 100644
--- a/scipy-stubs/sparse/_base.pyi
+++ b/scipy-stubs/sparse/_base.pyi
@@ -25,7 +25,6 @@ from ._typing import (
     Float,
     Index1D,
     Int,
-    Matrix,
     Scalar,
     Shape,
     ShapeCOO,
@@ -161,12 +160,12 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
     @overload  # dense array-like int
     def __init__(self: _spbase[np.int_], /, arg1: _ToSparseFromPy[opt.JustInt], *, maxprint: int | None = 50) -> None: ...
     @overload  # dense array-like float
-    def __init__(self: _spbase[np.float64], /, arg1: _ToSparseFromPy[opt.Just[float]], *, maxprint: int | None = 50) -> None: ...
+    def __init__(self: _spbase[np.float64], /, arg1: _ToSparseFromPy[opt.JustFloat], *, maxprint: int | None = 50) -> None: ...
     @overload  # dense array-like cfloat
     def __init__(
         self: _spbase[np.complex128],
         /,
-        arg1: _ToSparseFromPy[opt.Just[complex]],
+        arg1: _ToSparseFromPy[opt.JustComplex],
         *,
         maxprint: int | None = 50,
     ) -> None: ...
@@ -272,11 +271,11 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
     @overload  # sparse[-Complex], sparse[+Complex64]
     def __add__(self: _SpFromComplexT, other: _spbase[_ToComplex64], /) -> _SpFromComplexT: ...
     @overload  # spmatrix[-Int], array[+Int8]
-    def __add__(self: spmatrix[_FromInt], other: onp.ArrayND[_ToInt8], /) -> Matrix[_SCT_co]: ...
+    def __add__(self: spmatrix[_FromInt], other: onp.ArrayND[_ToInt8], /) -> onp.Matrix[_SCT_co]: ...
     @overload  # spmatrix[-Float], array[+Float32]
-    def __add__(self: spmatrix[_FromFloat], other: onp.ArrayND[_ToFloat32], /) -> Matrix[_SCT_co]: ...
+    def __add__(self: spmatrix[_FromFloat], other: onp.ArrayND[_ToFloat32], /) -> onp.Matrix[_SCT_co]: ...
     @overload  # spmatrix[-Complex], array[+Complex64]
-    def __add__(self: spmatrix[_FromComplex], other: onp.ArrayND[_ToComplex64], /) -> Matrix[_SCT_co]: ...
+    def __add__(self: spmatrix[_FromComplex], other: onp.ArrayND[_ToComplex64], /) -> onp.Matrix[_SCT_co]: ...
     @overload  # sparse[-Int], array[+Int8]
     def __add__(self: _spbase[_FromInt], other: onp.ArrayND[_ToInt8], /) -> onp.ArrayND[_SCT_co, _ShapeT_co]: ...
     @overload  # sparse[-Float], array[+Float32]
@@ -297,11 +296,11 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
     @overload  # sparse[-Complex], sparse[+Complex64]
     def __sub__(self: _SpFromComplexT, other: _spbase[_ToComplex64], /) -> _SpFromComplexT: ...
     @overload  # spmatrix[-Int], array[+Int8]
-    def __sub__(self: spmatrix[_FromInt], other: onp.ArrayND[_ToInt8], /) -> Matrix[_SCT_co]: ...
+    def __sub__(self: spmatrix[_FromInt], other: onp.ArrayND[_ToInt8], /) -> onp.Matrix[_SCT_co]: ...
     @overload  # spmatrix[-Float], array[+Float32]
-    def __sub__(self: spmatrix[_FromFloat], other: onp.ArrayND[_ToFloat32], /) -> Matrix[_SCT_co]: ...
+    def __sub__(self: spmatrix[_FromFloat], other: onp.ArrayND[_ToFloat32], /) -> onp.Matrix[_SCT_co]: ...
     @overload  # spmatrix[-Complex], array[+Complex64]
-    def __sub__(self: spmatrix[_FromComplex], other: onp.ArrayND[_ToComplex64], /) -> Matrix[_SCT_co]: ...
+    def __sub__(self: spmatrix[_FromComplex], other: onp.ArrayND[_ToComplex64], /) -> onp.Matrix[_SCT_co]: ...
     @overload  # sparse[-Int], array[+Int8]
     def __sub__(self: _spbase[_FromInt], other: onp.ArrayND[_ToInt8], /) -> onp.ArrayND[_SCT_co, _ShapeT_co]: ...
     @overload  # sparse[-Float], array[+Float32]
@@ -360,13 +359,13 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
     @overload  # spmatrix[+Bool], /, other: array-like ~Int
     def __mul__(self: spmatrix[_ToBool], /, other: _DualArrayLike2d[opt.JustInt, Int]) -> onp.Array2D[Int]: ...
     @overload  # spmatrix[+Int], /, other: scalar- or matrix-like ~Float
-    def __mul__(self: spmatrix[_ToInt], /, other: _DualMatrixLike[opt.Just[float], Float]) -> spmatrix[Float]: ...
+    def __mul__(self: spmatrix[_ToInt], /, other: _DualMatrixLike[opt.JustFloat, Float]) -> spmatrix[Float]: ...
     @overload  # spmatrix[+Int], /, other: array-like ~Float
-    def __mul__(self: spmatrix[_ToInt], /, other: _DualArrayLike2d[opt.Just[float], Float]) -> onp.Array2D[Float]: ...
+    def __mul__(self: spmatrix[_ToInt], /, other: _DualArrayLike2d[opt.JustFloat, Float]) -> onp.Array2D[Float]: ...
     @overload  # spmatrix[+Float], /, other: scalar- or matrix-like ~Complex
-    def __mul__(self: spmatrix[_ToFloat], /, other: _DualMatrixLike[opt.Just[complex], Complex]) -> spmatrix[Complex]: ...
+    def __mul__(self: spmatrix[_ToFloat], /, other: _DualMatrixLike[opt.JustComplex, Complex]) -> spmatrix[Complex]: ...
     @overload  # spmatrix[+Float], /, other: array-like ~Complex
-    def __mul__(self: spmatrix[_ToFloat], /, other: _DualArrayLike2d[opt.Just[complex], Complex]) -> onp.Array2D[Complex]: ...
+    def __mul__(self: spmatrix[_ToFloat], /, other: _DualArrayLike2d[opt.JustComplex, Complex]) -> onp.Array2D[Complex]: ...
     @overload  # Self[+Bool], /, other: -Int
     def __mul__(self: _spbase[_ToBool], /, other: _FromIntT) -> _spbase[_FromIntT, _ShapeT_co]: ...
     @overload  # Self[+Int], /, other: -Float
@@ -424,11 +423,11 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
     @overload  # spmatrix[+Int], other: scalar- or matrix-like ~Float
     def __matmul__(self: spmatrix[_ToInt], other: _spbase[Float], /) -> _SpMatrixOut[Float]: ...
     @overload  # spmatrix[+Int], other: array-like ~Float
-    def __matmul__(self: spmatrix[_ToInt], other: _DualArrayLike2d[opt.Just[float], Float], /) -> onp.Array2D[Float]: ...
+    def __matmul__(self: spmatrix[_ToInt], other: _DualArrayLike2d[opt.JustFloat, Float], /) -> onp.Array2D[Float]: ...
     @overload  # spmatrix[+Float], other: scalar- or matrix-like ~Complex
     def __matmul__(self: spmatrix[_ToFloat], other: _spbase[Complex], /) -> _SpMatrixOut[Complex]: ...
     @overload  # spmatrix[+Float], other: array-like ~Complex
-    def __matmul__(self: spmatrix[_ToFloat], other: _DualArrayLike2d[opt.Just[complex], Complex], /) -> onp.Array2D[Complex]: ...
+    def __matmul__(self: spmatrix[_ToFloat], other: _DualArrayLike2d[opt.JustComplex, Complex], /) -> onp.Array2D[Complex]: ...
     @overload  # catch-all
     def __matmul__(self, other: _DualArrayLike[complex, Scalar] | _spbase, /) -> _spbase[Any, Any] | onp.Array[Any, Any]: ...
     __rmatmul__ = __matmul__
@@ -449,11 +448,11 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
     @overload
     def __truediv__(self: spmatrix, rhs: onp.ArrayND[Scalar], /) -> coo_matrix[np.int_ | _FromFloat]: ...
     @overload
-    def __truediv__(self: spmatrix, rhs: _spbase[_FromComplex], /) -> Matrix[_FromComplex]: ...
+    def __truediv__(self: spmatrix, rhs: _spbase[_FromComplex], /) -> onp.Matrix[_FromComplex]: ...
     @overload
-    def __truediv__(self: spmatrix, rhs: _spbase[_FromFloat], /) -> Matrix[_FromFloat]: ...
+    def __truediv__(self: spmatrix, rhs: _spbase[_FromFloat], /) -> onp.Matrix[_FromFloat]: ...
     @overload
-    def __truediv__(self: spmatrix, rhs: _spbase, /) -> Matrix[np.int_ | _FromFloat]: ...
+    def __truediv__(self: spmatrix, rhs: _spbase, /) -> onp.Matrix[np.int_ | _FromFloat]: ...
     @overload
     def __truediv__(self, rhs: _spbase[_FromComplex], /) -> onp.ArrayND[_FromComplex, _ShapeT_co]: ...
     @overload
@@ -477,9 +476,9 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
     @overload  # spmatrix[+Bool], /, other: scalar-like ~Int
     def dot(self: spmatrix[_ToBool], /, other: opt.JustInt) -> spmatrix[Int]: ...
     @overload  # spmatrix[+Int], /, other: scalar-like ~Float
-    def dot(self: spmatrix[_ToInt], /, other: opt.Just[float]) -> spmatrix[Float]: ...
+    def dot(self: spmatrix[_ToInt], /, other: opt.JustFloat) -> spmatrix[Float]: ...
     @overload  # spmatrix[+Float], /, other: scalar-like ~Complex
-    def dot(self: spmatrix[_ToFloat], /, other: opt.Just[complex]) -> spmatrix[Complex]: ...
+    def dot(self: spmatrix[_ToFloat], /, other: opt.JustComplex) -> spmatrix[Complex]: ...
     @overload  # sparray[-Bool], /, other: sparse +Bool
     def dot(self: _SpArray, /, other: _spbase[_ToBool | _SCT_co]) -> _SpArrayOut[_SCT_co]: ...
     @overload  # sparray[-Bool], /, other: array-like +Bool
@@ -521,11 +520,11 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
     @overload  # spmatrix[+Int], /, other: scalar- or matrix-like ~Float
     def dot(self: spmatrix[_ToInt], /, other: _spbase[Float]) -> _SpMatrixOut[Float]: ...
     @overload  # spmatrix[+Int], /, other: array-like ~Float
-    def dot(self: spmatrix[_ToInt], /, other: _DualArrayLike2d[opt.Just[float], Float]) -> onp.Array2D[Float]: ...
+    def dot(self: spmatrix[_ToInt], /, other: _DualArrayLike2d[opt.JustFloat, Float]) -> onp.Array2D[Float]: ...
     @overload  # spmatrix[+Float], /, other: scalar- or matrix-like ~Complex
     def dot(self: spmatrix[_ToFloat], /, other: _spbase[Complex]) -> _SpMatrixOut[Complex]: ...
     @overload  # spmatrix[+Float], /, other: array-like ~Complex
-    def dot(self: spmatrix[_ToFloat], /, other: _DualArrayLike2d[opt.Just[complex], Complex]) -> onp.Array2D[Complex]: ...
+    def dot(self: spmatrix[_ToFloat], /, other: _DualArrayLike2d[opt.JustComplex, Complex]) -> onp.Array2D[Complex]: ...
     @overload  # catch-all
     def dot(self, /, other: _DualArrayLike[complex, Scalar] | _spbase) -> _spbase[Any, Any] | onp.Array[Any, Any]: ...
 
@@ -551,7 +550,7 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
     @overload  # sparray, dense
     def maximum(self: sparray, /, other: onp.ArrayND[_SCT]) -> onp.ArrayND[_SCT_co | _SCT, _ShapeT_co]: ...
     @overload  # spmatrix, dense
-    def maximum(self: spmatrix, /, other: onp.ArrayND[_SCT]) -> Matrix[_SCT_co | _SCT]: ...
+    def maximum(self: spmatrix, /, other: onp.ArrayND[_SCT]) -> onp.Matrix[_SCT_co | _SCT]: ...
     @overload  # sparray, unknown scalar
     def maximum(self: sparray, /, other: onp.ToComplex) -> _SpArrayOut[Any]: ...
     @overload  # spmatrix, unknown scalar
@@ -579,7 +578,7 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
     @overload  # sparray, dense
     def minimum(self: sparray, /, other: onp.ArrayND[_SCT]) -> onp.ArrayND[_SCT_co | _SCT, _ShapeT_co]: ...
     @overload  # spmatrix, dense
-    def minimum(self: spmatrix, /, other: onp.ArrayND[_SCT]) -> Matrix[_SCT_co | _SCT]: ...
+    def minimum(self: spmatrix, /, other: onp.ArrayND[_SCT]) -> onp.Matrix[_SCT_co | _SCT]: ...
     @overload  # sparray, unknown scalar
     def minimum(self: sparray, /, other: onp.ToComplex) -> _SpArrayOut[Any]: ...
     @overload  # spmatrix, unknown scalar
@@ -622,7 +621,7 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
     @overload  # sparray[-Bool, 2d], axis: index
     def sum(self: _SpArray2D[np.bool_], /, axis: op.CanIndex, dtype: None = None, out: None = None) -> onp.Array1D[np.int_]: ...
     @overload  # spmatrix[-Bool], axis: index
-    def sum(self: spmatrix[np.bool_], /, axis: op.CanIndex, dtype: None = None, out: None = None) -> Matrix[np.int_]: ...
+    def sum(self: spmatrix[np.bool_], /, axis: op.CanIndex, dtype: None = None, out: None = None) -> onp.Matrix[np.int_]: ...
     @overload  # Self[-Int]
     def sum(self: _spbase[_FromIntT], /, axis: None = None, dtype: None = None, out: None = None) -> _FromIntT: ...
     @overload  # sparray[-Int, 1d], axis: index
@@ -642,9 +641,9 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
         out: None = None,
     ) -> onp.Array1D[_FromIntT]: ...
     @overload  # spmatrix[-Int], axis: index
-    def sum(self: spmatrix[_FromIntT], /, axis: op.CanIndex, dtype: None = None, out: None = None) -> Matrix[_FromIntT]: ...
+    def sum(self: spmatrix[_FromIntT], /, axis: op.CanIndex, dtype: None = None, out: None = None) -> onp.Matrix[_FromIntT]: ...
     @overload  # spmatrix, axis: index, dtype: <unknown>
-    def sum(self: spmatrix, /, axis: op.CanIndex, dtype: npt.DTypeLike, out: None = None) -> Matrix[Any]: ...
+    def sum(self: spmatrix, /, axis: op.CanIndex, dtype: npt.DTypeLike, out: None = None) -> onp.Matrix[Any]: ...
     @overload  # dtype: <unknown>  (keyword)
     def sum(self, /, axis: op.CanIndex | None = None, *, dtype: npt.DTypeLike, out: None = None) -> Any: ...  # noqa: ANN401
 
@@ -680,7 +679,7 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
         axis: op.CanIndex,
         dtype: None = None,
         out: None = None,
-    ) -> Matrix[np.float64]: ...
+    ) -> onp.Matrix[np.float64]: ...
     @overload  # sparray[-Float, 1d]
     def mean(
         self: _SpArray1D[_FromFloatT],
@@ -698,9 +697,11 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
         out: None = None,
     ) -> onp.Array1D[_FromFloatT]: ...
     @overload  # spmatrix[-Float], axis: index
-    def mean(self: spmatrix[_FromFloatT], /, axis: op.CanIndex, dtype: None = None, out: None = None) -> Matrix[_FromFloatT]: ...
+    def mean(
+        self: spmatrix[_FromFloatT], /, axis: op.CanIndex, dtype: None = None, out: None = None
+    ) -> onp.Matrix[_FromFloatT]: ...
     @overload  # spmatrix, axis: index, dtype: <unknown>
-    def mean(self: spmatrix, /, axis: op.CanIndex, dtype: npt.DTypeLike, out: None = None) -> Matrix[Any]: ...
+    def mean(self: spmatrix, /, axis: op.CanIndex, dtype: npt.DTypeLike, out: None = None) -> onp.Matrix[Any]: ...
     @overload  # dtype: <unknown>  (keyword)
     def mean(self, /, axis: op.CanIndex | None = None, *, dtype: npt.DTypeLike, out: None = None) -> Any: ...  # noqa: ANN401
 
@@ -925,11 +926,11 @@ class _spbase(Generic[_SCT_co, _ShapeT_co]):
 
     #
     @overload  # self: spmatrix, out: None
-    def todense(self: spmatrix, /, order: OrderCF | None = None, out: None = None) -> Matrix[_SCT_co]: ...
+    def todense(self: spmatrix, /, order: OrderCF | None = None, out: None = None) -> onp.Matrix[_SCT_co]: ...
     @overload  # self: spmatrix, out: array (positional)
-    def todense(self: spmatrix, /, order: OrderCF | None, out: onp.ArrayND[_SCT]) -> Matrix[_SCT]: ...
+    def todense(self: spmatrix, /, order: OrderCF | None, out: onp.ArrayND[_SCT]) -> onp.Matrix[_SCT]: ...
     @overload  # self: spmatrix, out: array (keyword)
-    def todense(self: spmatrix, /, order: OrderCF | None = None, *, out: onp.ArrayND[_SCT]) -> Matrix[_SCT]: ...
+    def todense(self: spmatrix, /, order: OrderCF | None = None, *, out: onp.ArrayND[_SCT]) -> onp.Matrix[_SCT]: ...
     @overload  # out: None
     def todense(self, /, order: OrderCF | None = None, out: None = None) -> onp.Array[_ShapeT_co, _SCT_co]: ...
     @overload  # out: array (positional)
diff --git a/scipy-stubs/sparse/_bsr.pyi b/scipy-stubs/sparse/_bsr.pyi
index 5f1862e8..3fee02b7 100644
--- a/scipy-stubs/sparse/_bsr.pyi
+++ b/scipy-stubs/sparse/_bsr.pyi
@@ -92,9 +92,9 @@ class _bsr_base(_cs_matrix[_SCT, tuple[int, int]], _minmax_mixin[_SCT, tuple[int
     def __init__(
         self: _bsr_base[np.float64],
         /,
-        arg1: _ToMatrixPy[opt.Just[float]],
+        arg1: _ToMatrixPy[opt.JustFloat],
         shape: ToShape2d | None = None,
-        dtype: type[opt.Just[float]] | onp.AnyFloat64DType | None = None,
+        dtype: type[opt.JustFloat] | onp.AnyFloat64DType | None = None,
         copy: bool = False,
         blocksize: tuple[int, int] | None = None,
         *,
@@ -104,9 +104,9 @@ class _bsr_base(_cs_matrix[_SCT, tuple[int, int]], _minmax_mixin[_SCT, tuple[int
     def __init__(
         self: _bsr_base[np.complex128],
         /,
-        arg1: _ToMatrixPy[opt.Just[complex]],
+        arg1: _ToMatrixPy[opt.JustComplex],
         shape: ToShape2d | None = None,
-        dtype: type[opt.Just[complex]] | onp.AnyComplex128DType | None = None,
+        dtype: type[opt.JustComplex] | onp.AnyComplex128DType | None = None,
         copy: bool = False,
         blocksize: tuple[int, int] | None = None,
         *,
diff --git a/scipy-stubs/sparse/_compressed.pyi b/scipy-stubs/sparse/_compressed.pyi
index 81d521bf..441aea27 100644
--- a/scipy-stubs/sparse/_compressed.pyi
+++ b/scipy-stubs/sparse/_compressed.pyi
@@ -184,7 +184,7 @@ class _cs_matrix(
     def __init__(
         self: _cs_matrix[np.float64, tuple[int]],
         /,
-        arg1: Sequence[opt.Just[float]],
+        arg1: Sequence[opt.JustFloat],
         shape: ToShape1dNd | None = None,
         dtype: ToDTypeFloat | None = None,
         copy: bool = False,
@@ -195,7 +195,7 @@ class _cs_matrix(
     def __init__(
         self: _cs_matrix[np.float64, tuple[int, int]],
         /,
-        arg1: Sequence[Sequence[opt.Just[float]]],
+        arg1: Sequence[Sequence[opt.JustFloat]],
         shape: ToShape1dNd | None = None,
         dtype: ToDTypeFloat | None = None,
         copy: bool = False,
@@ -206,7 +206,7 @@ class _cs_matrix(
     def __init__(
         self: _cs_matrix[np.complex128, tuple[int]],
         /,
-        arg1: Sequence[opt.Just[complex]],
+        arg1: Sequence[opt.JustComplex],
         shape: ToShape1dNd | None = None,
         dtype: ToDTypeComplex | None = None,
         copy: bool = False,
@@ -217,7 +217,7 @@ class _cs_matrix(
     def __init__(
         self: _cs_matrix[np.complex128, tuple[int, int]],
         /,
-        arg1: Sequence[Sequence[opt.Just[complex]]],
+        arg1: Sequence[Sequence[opt.JustComplex]],
         shape: ToShape1dNd | None = None,
         dtype: ToDTypeComplex | None = None,
         copy: bool = False,
diff --git a/scipy-stubs/sparse/_construct.pyi b/scipy-stubs/sparse/_construct.pyi
index 0dc42c57..ac047896 100644
--- a/scipy-stubs/sparse/_construct.pyi
+++ b/scipy-stubs/sparse/_construct.pyi
@@ -139,9 +139,9 @@ class _DataSampler(Protocol):
 
 ###
 
-#
+# NOTE: The `overload-overlap` mypy errors are false positives.
 @overload  # diagonals: <known>, dtype: None = ..., format: {"dia", None} = ...
-def diags_array(
+def diags_array(  # type: ignore[overload-overlap]
     diagonals: _ToArray1D[_SCT] | _ToArray2D[_SCT],
     /,
     *,
@@ -151,7 +151,7 @@ def diags_array(
     dtype: None = None,
 ) -> dia_array[_SCT]: ...
 @overload  # diagonals: <known>, dtype: None = ..., format: <otherwise>
-def diags_array(
+def diags_array(  # type: ignore[overload-overlap]
     diagonals: _ToArray1D[_SCT] | _ToArray2D[_SCT],
     /,
     *,
@@ -652,7 +652,7 @@ def kron(
     format: SPFormat | None = None,
 ) -> _SpArray2D[_SCT] | _SpMatrix[_SCT]: ...
 
-#
+# NOTE: The `overload-overlap` mypy errors are false positives.
 @overload  # A: spmatrix or 2d array-like, B: spmatrix or 2d array-like, format: {"csr", None} = ...
 def kronsum(A: _ToSpMatrix[_SCT1], B: _ToSpMatrix[_SCT2], format: _FmtCSR | None = None) -> csr_matrix[_SCT1 | _SCT2]: ...
 @overload  # A: spmatrix or 2d array-like, B: spmatrix or 2d array-like, format: <otherwise>
@@ -702,9 +702,9 @@ def hstack(blocks: Seq[_SpArray], format: _FmtOut | None = None, *, dtype: npt.D
 @overload  # sparray, format: <non-default>, dtype: <unknown>
 def hstack(blocks: Seq[_SpArray], format: _FmtNonOut, dtype: npt.DTypeLike) -> _SpArrayNonOut: ...
 @overload  # spmatrix, format: <default>, dtype: None
-def hstack(blocks: Seq[spmatrix[_SCT]], format: _FmtOut | None = None, dtype: None = None) -> _SpMatrixOut[_SCT]: ...
+def hstack(blocks: Seq[spmatrix[_SCT]], format: _FmtOut | None = None, dtype: None = None) -> _SpMatrixOut[_SCT]: ...  # type: ignore[overload-overlap]
 @overload  # spmatrix, format: <non-default>, dtype: None
-def hstack(blocks: Seq[spmatrix[_SCT]], format: _FmtNonOut, dtype: None = None) -> _SpMatrixNonOut[_SCT]: ...
+def hstack(blocks: Seq[spmatrix[_SCT]], format: _FmtNonOut, dtype: None = None) -> _SpMatrixNonOut[_SCT]: ...  # type: ignore[overload-overlap]
 @overload  # spmatrix, format: <default>, dtype: <int>
 def hstack(blocks: Seq[spmatrix], format: _FmtOut | None = None, *, dtype: ToDTypeBool) -> _SpMatrixOut[np.bool_]: ...
 @overload  # spmatrix, format: <non-default>, dtype: <int>
@@ -760,9 +760,9 @@ def vstack(blocks: Seq[_SpArray], format: _FmtOut | None = None, *, dtype: npt.D
 @overload  # sparray, format: <non-default>, dtype: <unknown>
 def vstack(blocks: Seq[_SpArray], format: _FmtNonOut, dtype: npt.DTypeLike) -> _SpArrayNonOut: ...
 @overload  # spmatrix, format: <default>, dtype: None
-def vstack(blocks: Seq[spmatrix[_SCT]], format: _FmtOut | None = None, dtype: None = None) -> _SpMatrixOut[_SCT]: ...
+def vstack(blocks: Seq[spmatrix[_SCT]], format: _FmtOut | None = None, dtype: None = None) -> _SpMatrixOut[_SCT]: ...  # type: ignore[overload-overlap]
 @overload  # spmatrix, format: <non-default>, dtype: None
-def vstack(blocks: Seq[spmatrix[_SCT]], format: _FmtNonOut, dtype: None = None) -> _SpMatrixNonOut[_SCT]: ...
+def vstack(blocks: Seq[spmatrix[_SCT]], format: _FmtNonOut, dtype: None = None) -> _SpMatrixNonOut[_SCT]: ...  # type: ignore[overload-overlap]
 @overload  # spmatrix, format: <default>, dtype: <int>
 def vstack(blocks: Seq[spmatrix], format: _FmtOut | None = None, *, dtype: ToDTypeBool) -> _SpMatrixOut[np.bool_]: ...
 @overload  # spmatrix, format: <non-default>, dtype: <int>
diff --git a/scipy-stubs/sparse/_coo.pyi b/scipy-stubs/sparse/_coo.pyi
index d23852e5..8b85697b 100644
--- a/scipy-stubs/sparse/_coo.pyi
+++ b/scipy-stubs/sparse/_coo.pyi
@@ -183,7 +183,7 @@ class _coo_base(_data_matrix[_SCT, _ShapeT_co], _minmax_mixin[_SCT, _ShapeT_co],
     def __init__(
         self: _coo_base[np.float64, tuple[int]],
         /,
-        arg1: Sequence[opt.Just[float]],
+        arg1: Sequence[opt.JustFloat],
         shape: ToShape1d | None = None,
         dtype: ToDTypeFloat | None = None,
         copy: bool = False,
@@ -194,7 +194,7 @@ class _coo_base(_data_matrix[_SCT, _ShapeT_co], _minmax_mixin[_SCT, _ShapeT_co],
     def __init__(
         self: _coo_base[np.float64, tuple[int, int]],
         /,
-        arg1: Sequence[Sequence[opt.Just[float]]],
+        arg1: Sequence[Sequence[opt.JustFloat]],
         shape: ToShape2d | None = None,
         dtype: ToDTypeFloat | None = None,
         copy: bool = False,
@@ -205,7 +205,7 @@ class _coo_base(_data_matrix[_SCT, _ShapeT_co], _minmax_mixin[_SCT, _ShapeT_co],
     def __init__(
         self: _coo_base[np.complex128, tuple[int]],
         /,
-        arg1: Sequence[opt.Just[complex]],
+        arg1: Sequence[opt.JustComplex],
         shape: ToShape1d | None = None,
         dtype: ToDTypeComplex | None = None,
         copy: bool = False,
@@ -216,7 +216,7 @@ class _coo_base(_data_matrix[_SCT, _ShapeT_co], _minmax_mixin[_SCT, _ShapeT_co],
     def __init__(
         self: _coo_base[np.complex128, tuple[int, int]],
         /,
-        arg1: Sequence[Sequence[opt.Just[complex]]],
+        arg1: Sequence[Sequence[opt.JustComplex]],
         shape: ToShape2d | None = None,
         dtype: ToDTypeComplex | None = None,
         copy: bool = False,
diff --git a/scipy-stubs/sparse/_data.pyi b/scipy-stubs/sparse/_data.pyi
index 4983391c..4a348f10 100644
--- a/scipy-stubs/sparse/_data.pyi
+++ b/scipy-stubs/sparse/_data.pyi
@@ -9,7 +9,7 @@ from ._base import _spbase, sparray
 from ._coo import coo_array, coo_matrix
 from ._matrix import spmatrix
 from ._sputils import _ScalarLike
-from ._typing import Matrix, Scalar, ShapeCOO, ToShape1d, ToShape2d
+from ._typing import Scalar, ShapeCOO, ToShape1d, ToShape2d
 
 __all__: list[str] = []
 
@@ -250,7 +250,7 @@ class _minmax_mixin(Generic[_SCT_co, _ShapeT_co]):
     @overload  # sparray, axis: int, out: None = ...
     def argmax(self: sparray, /, axis: onp.ToInt, out: None = None, *, explicit: bool = False) -> onp.Array1D[np.intp]: ...  # type: ignore[misc]
     @overload  # spmatrix, axis: int, out: None = ...
-    def argmax(self: spmatrix, /, axis: onp.ToInt, out: None = None, *, explicit: bool = False) -> Matrix[np.intp]: ...  # type: ignore[misc]
+    def argmax(self: spmatrix, /, axis: onp.ToInt, out: None = None, *, explicit: bool = False) -> onp.Matrix[np.intp]: ...  # type: ignore[misc]
 
     #
     @overload  # axis: None = ..., out: None = ...
@@ -267,4 +267,4 @@ class _minmax_mixin(Generic[_SCT_co, _ShapeT_co]):
     @overload  # sparray, axis: int, out: None = ...
     def argmin(self: sparray, /, axis: onp.ToInt, out: None = None, *, explicit: bool = False) -> onp.Array1D[np.intp]: ...  # type: ignore[misc]
     @overload  # spmatrix, axis: int, out: None = ...
-    def argmin(self: spmatrix, /, axis: onp.ToInt, out: None = None, *, explicit: bool = False) -> Matrix[np.intp]: ...  # type: ignore[misc]
+    def argmin(self: spmatrix, /, axis: onp.ToInt, out: None = None, *, explicit: bool = False) -> onp.Matrix[np.intp]: ...  # type: ignore[misc]
diff --git a/scipy-stubs/sparse/_dia.pyi b/scipy-stubs/sparse/_dia.pyi
index 196cd702..ae0ef391 100644
--- a/scipy-stubs/sparse/_dia.pyi
+++ b/scipy-stubs/sparse/_dia.pyi
@@ -85,9 +85,9 @@ class _dia_base(_data_matrix[_SCT, tuple[int, int]], Generic[_SCT]):
     def __init__(
         self: _dia_base[np.float64],
         /,
-        arg1: _ToMatrixPy[opt.Just[float]],
+        arg1: _ToMatrixPy[opt.JustFloat],
         shape: ToShape2d | None = None,
-        dtype: type[opt.Just[float]] | onp.AnyFloat64DType | None = None,
+        dtype: type[opt.JustFloat] | onp.AnyFloat64DType | None = None,
         copy: bool = False,
         *,
         maxprint: int | None = None,
@@ -96,9 +96,9 @@ class _dia_base(_data_matrix[_SCT, tuple[int, int]], Generic[_SCT]):
     def __init__(
         self: _dia_base[np.complex128],
         /,
-        arg1: _ToMatrixPy[opt.Just[complex]],
+        arg1: _ToMatrixPy[opt.JustComplex],
         shape: ToShape2d | None = None,
-        dtype: type[opt.Just[complex]] | onp.AnyComplex128DType | None = None,
+        dtype: type[opt.JustComplex] | onp.AnyComplex128DType | None = None,
         copy: bool = False,
         *,
         maxprint: int | None = None,
diff --git a/scipy-stubs/sparse/_dok.pyi b/scipy-stubs/sparse/_dok.pyi
index 6b31d0ff..70866560 100644
--- a/scipy-stubs/sparse/_dok.pyi
+++ b/scipy-stubs/sparse/_dok.pyi
@@ -89,9 +89,9 @@ class _dok_base(_spbase[_SCT, _ShapeT_co], IndexMixin[_SCT, _ShapeT_co], dict[Sh
     def __init__(
         self: _dok_base[np.float64],
         /,
-        arg1: _ToMatrixPy[opt.Just[float]],
+        arg1: _ToMatrixPy[opt.JustFloat],
         shape: ToShape1dNd | None = None,
-        dtype: type[opt.Just[float]] | onp.AnyFloat64DType | None = None,
+        dtype: type[opt.JustFloat] | onp.AnyFloat64DType | None = None,
         copy: bool = False,
         *,
         maxprint: int | None = None,
@@ -100,9 +100,9 @@ class _dok_base(_spbase[_SCT, _ShapeT_co], IndexMixin[_SCT, _ShapeT_co], dict[Sh
     def __init__(
         self: _dok_base[np.complex128],
         /,
-        arg1: _ToMatrixPy[opt.Just[complex]],
+        arg1: _ToMatrixPy[opt.JustComplex],
         shape: ToShape1dNd | None = None,
-        dtype: type[opt.Just[complex]] | onp.AnyComplex128DType | None = None,
+        dtype: type[opt.JustComplex] | onp.AnyComplex128DType | None = None,
         copy: bool = False,
         *,
         maxprint: int | None = None,
diff --git a/scipy-stubs/sparse/_lil.pyi b/scipy-stubs/sparse/_lil.pyi
index 46d37f9b..ef73497e 100644
--- a/scipy-stubs/sparse/_lil.pyi
+++ b/scipy-stubs/sparse/_lil.pyi
@@ -87,7 +87,7 @@ class _lil_base(_spbase[_SCT, tuple[int, int]], IndexMixin[_SCT, tuple[int, int]
     def __init__(
         self: _lil_base[np.float64],
         /,
-        arg1: _ToMatrixPy[opt.Just[float]],
+        arg1: _ToMatrixPy[opt.JustFloat],
         shape: ToShape2d | None = None,
         dtype: ToDTypeFloat | None = None,
         copy: bool = False,
@@ -98,7 +98,7 @@ class _lil_base(_spbase[_SCT, tuple[int, int]], IndexMixin[_SCT, tuple[int, int]
     def __init__(
         self: _lil_base[np.complex128],
         /,
-        arg1: _ToMatrixPy[opt.Just[complex]],
+        arg1: _ToMatrixPy[opt.JustComplex],
         shape: ToShape2d | None = None,
         dtype: ToDTypeComplex | None = None,
         copy: bool = False,
diff --git a/scipy-stubs/sparse/_matrix.pyi b/scipy-stubs/sparse/_matrix.pyi
index 6ef57c20..38ff3a21 100644
--- a/scipy-stubs/sparse/_matrix.pyi
+++ b/scipy-stubs/sparse/_matrix.pyi
@@ -20,7 +20,7 @@ from ._csr import csr_matrix
 from ._dia import dia_matrix
 from ._dok import dok_matrix
 from ._lil import lil_matrix
-from ._typing import Complex, Float, Int, Matrix, Scalar, SPFormat, ToShape2d
+from ._typing import Complex, Float, Int, Scalar, SPFormat, ToShape2d
 
 _T = TypeVar("_T")
 _SCT = TypeVar("_SCT", bound=Scalar)
@@ -104,13 +104,13 @@ class spmatrix(Generic[_SCT_co]):
     @overload  # spmatrix[+Bool], other: array-like ~Int
     def __mul__(self: spmatrix[_ToBool], other: _DualArrayLike[opt.JustInt, Int], /) -> onp.Array2D[Int]: ...
     @overload  # spmatrix[+Int], other: scalar- or matrix-like ~Float
-    def __mul__(self: spmatrix[_ToInt], other: _DualMatrixLike[opt.Just[float], Float], /) -> spmatrix[Float]: ...
+    def __mul__(self: spmatrix[_ToInt], other: _DualMatrixLike[opt.JustFloat, Float], /) -> spmatrix[Float]: ...
     @overload  # spmatrix[+Int], other: array-like ~Float
-    def __mul__(self: spmatrix[_ToInt], other: _DualArrayLike[opt.Just[float], Float], /) -> onp.Array2D[Float]: ...
+    def __mul__(self: spmatrix[_ToInt], other: _DualArrayLike[opt.JustFloat, Float], /) -> onp.Array2D[Float]: ...
     @overload  # spmatrix[+Float], other: scalar- or matrix-like ~Complex
-    def __mul__(self: spmatrix[_ToFloat], other: _DualMatrixLike[opt.Just[complex], Complex], /) -> spmatrix[Complex]: ...
+    def __mul__(self: spmatrix[_ToFloat], other: _DualMatrixLike[opt.JustComplex, Complex], /) -> spmatrix[Complex]: ...
     @overload  # spmatrix[+Float], other: array-like ~Complex
-    def __mul__(self: spmatrix[_ToFloat], other: _DualArrayLike[opt.Just[complex], Complex], /) -> onp.Array2D[Complex]: ...
+    def __mul__(self: spmatrix[_ToFloat], other: _DualArrayLike[opt.JustComplex, Complex], /) -> onp.Array2D[Complex]: ...
     @overload  # catch-all
     def __mul__(self, other: _DualArrayLike[complex, Scalar] | _spbase, /) -> _spbase[Any, Any] | onp.Array[Any, Any]: ...
     __rmul__ = __mul__
@@ -137,8 +137,8 @@ class spmatrix(Generic[_SCT_co]):
 
     #
     @overload
-    def todense(self, /, order: OrderCF | None = None, out: None = None) -> Matrix[_SCT_co]: ...
+    def todense(self, /, order: OrderCF | None = None, out: None = None) -> onp.Matrix[_SCT_co]: ...
     @overload
-    def todense(self, /, order: OrderCF | None, out: onp.ArrayND[_SCT]) -> Matrix[_SCT]: ...
+    def todense(self, /, order: OrderCF | None, out: onp.ArrayND[_SCT]) -> onp.Matrix[_SCT]: ...
     @overload
-    def todense(self, /, order: OrderCF | None = None, *, out: onp.ArrayND[_SCT]) -> Matrix[_SCT]: ...
+    def todense(self, /, order: OrderCF | None = None, *, out: onp.ArrayND[_SCT]) -> onp.Matrix[_SCT]: ...
diff --git a/scipy-stubs/sparse/_sputils.pyi b/scipy-stubs/sparse/_sputils.pyi
index 2a5e0306..d10b0498 100644
--- a/scipy-stubs/sparse/_sputils.pyi
+++ b/scipy-stubs/sparse/_sputils.pyi
@@ -19,7 +19,7 @@ from scipy.sparse import (
     dia_array,
     dia_matrix,
 )
-from scipy.sparse._typing import Matrix, ToDType
+from scipy.sparse._typing import ToDType
 
 __all__ = [
     "broadcast_shapes",
@@ -143,10 +143,10 @@ def matrix(
     subok: bool = False,
     ndmin: Literal[0, 1, 2] = 0,
     like: onp.CanArrayFunction | None = None,
-) -> Matrix[_SCT]: ...
+) -> onp.Matrix[_SCT]: ...
 
 #
-def asmatrix(data: _ToArray2D[_SCT], dtype: ToDType[_SCT] | type | str | None = None) -> Matrix[_SCT]: ...
+def asmatrix(data: _ToArray2D[_SCT], dtype: ToDType[_SCT] | type | str | None = None) -> onp.Matrix[_SCT]: ...
 
 #
 @overload  # BSR/CSC/CSR, dtype: <default>
diff --git a/scipy-stubs/sparse/_typing.pyi b/scipy-stubs/sparse/_typing.pyi
index 191a85da..df818378 100644
--- a/scipy-stubs/sparse/_typing.pyi
+++ b/scipy-stubs/sparse/_typing.pyi
@@ -12,7 +12,6 @@ __all__ = (
     "Float",
     "Index1D",
     "Int",
-    "Matrix",
     "SPFormat",
     "Scalar",
     "Shape",
@@ -49,21 +48,16 @@ Complex: TypeAlias = np.complexfloating[Any, Any]
 # NOTE: The `TypeAliasType` is used to avoid long error messages.
 Scalar = TypeAliasType("Scalar", np.bool_ | Int | Float | Complex)
 
-Matrix: TypeAlias = np.matrix[tuple[int, int], np.dtype[_SCT]]
 Index1D: TypeAlias = onp.Array1D[np.int32 | np.int64]
 
 SPFormat: TypeAlias = Literal["bsr", "coo", "csc", "csr", "dia", "dok", "lil"]
 
 _SCT = TypeVar("_SCT", bound=np.generic)
 ToDType: TypeAlias = type[_SCT] | np.dtype[_SCT] | onp.HasDType[np.dtype[_SCT]]
-ToDTypeBool: TypeAlias = type[bool] | ToDType[np.bool_] | Literal["bool", "bool_", "?", "b1"]
-ToDTypeInt: TypeAlias = type[opt.JustInt] | ToDType[np.int_] | Literal["int", "int_", "n"]
-ToDTypeFloat: TypeAlias = type[opt.Just[float]] | ToDType[np.float64] | Literal["float", "float64", "double", "f8", "d"]
-ToDTypeComplex: TypeAlias = (
-    type[opt.Just[complex]]
-    | ToDType[np.complex128]
-    | Literal["complex", "complex128", "cdouble", "c16", "D"]
-)  # fmt: skip
+ToDTypeBool: TypeAlias = onp.AnyBoolDType
+ToDTypeInt: TypeAlias = type[opt.JustInt] | onp.AnyIntDType  # see https://github.com/jorenham/optype/issues/235
+ToDTypeFloat: TypeAlias = onp.AnyFloat64DType
+ToDTypeComplex: TypeAlias = onp.AnyComplex128DType
 
 ToShape1d: TypeAlias = tuple[op.CanIndex]  # ndim == 1
 ToShape2d: TypeAlias = tuple[op.CanIndex, op.CanIndex]  # ndim == 2
diff --git a/scipy-stubs/sparse/csgraph/_tools.pyi b/scipy-stubs/sparse/csgraph/_tools.pyi
index fa04c079..e7046d91 100644
--- a/scipy-stubs/sparse/csgraph/_tools.pyi
+++ b/scipy-stubs/sparse/csgraph/_tools.pyi
@@ -10,10 +10,11 @@ from scipy.sparse._typing import Float, Int
 _Real: TypeAlias = Int | Float
 _RealT = TypeVar("_RealT", bound=_Real, default=_Real)
 
-_MaskedArray2D: TypeAlias = np.ma.MaskedArray[tuple[int, int], np.dtype[_RealT]]
 _SparseGraph: TypeAlias = (
-    csr_array[_RealT] | csr_matrix[_RealT] | csc_array[_RealT] | csc_matrix[_RealT] | lil_array[_RealT] | lil_matrix[_RealT]
-)
+    csr_array[_RealT] | csr_matrix[_RealT]
+    | csc_array[_RealT] | csc_matrix[_RealT]
+    | lil_array[_RealT] | lil_matrix[_RealT]
+)  # fmt: skip
 
 _ToGraph: TypeAlias = onp.ToFloat2D | _spbase[_Real, tuple[int, int]]
 
@@ -22,7 +23,7 @@ _ToGraph: TypeAlias = onp.ToFloat2D | _spbase[_Real, tuple[int, int]]
 DTYPE: Final[type[np.float64]] = ...
 ITYPE: Final[type[np.int32]] = ...
 
-def csgraph_from_masked(graph: _MaskedArray2D) -> csr_array[np.float64 | np.int32]: ...
+def csgraph_from_masked(graph: onp.MArray2D[_Real]) -> csr_array[np.float64 | np.int32]: ...
 
 #
 def csgraph_masked_from_dense(
@@ -31,7 +32,7 @@ def csgraph_masked_from_dense(
     nan_null: bool = True,
     infinity_null: bool = True,
     copy: bool = True,
-) -> _MaskedArray2D[np.float64 | np.int32]: ...
+) -> onp.MArray2D[np.float64 | np.int32]: ...
 
 #
 def csgraph_from_dense(
@@ -43,7 +44,7 @@ def csgraph_from_dense(
 
 #
 def csgraph_to_dense(csgraph: _SparseGraph[_Real], null_value: float | None = 0) -> onp.Array2D[np.float64 | np.int32]: ...
-def csgraph_to_masked(csgraph: _SparseGraph[_Real]) -> _MaskedArray2D[np.float64 | np.int32]: ...
+def csgraph_to_masked(csgraph: _SparseGraph[_Real]) -> onp.MArray2D[np.float64 | np.int32]: ...
 
 #
 def reconstruct_path(
diff --git a/scipy-stubs/sparse/linalg/_interface.pyi b/scipy-stubs/sparse/linalg/_interface.pyi
index d0dae9cf..1435b050 100644
--- a/scipy-stubs/sparse/linalg/_interface.pyi
+++ b/scipy-stubs/sparse/linalg/_interface.pyi
@@ -10,7 +10,7 @@ import numpy as np
 import optype as op
 import optype.numpy as onp
 from scipy.sparse._base import _spbase
-from scipy.sparse._typing import Complex, Float, Int, Matrix, Scalar, ToDTypeComplex, ToDTypeFloat, ToDTypeInt
+from scipy.sparse._typing import Complex, Float, Int, Scalar, ToDTypeComplex, ToDTypeFloat, ToDTypeInt
 
 __all__ = ["LinearOperator", "aslinearoperator"]
 
@@ -70,9 +70,9 @@ class LinearOperator(Generic[_SCT_co]):
     @overload  # float array 1d
     def matvec(self, /, x: onp.ToFloatStrict1D) -> onp.Array1D[_SCT_co]: ...
     @overload  # float matrix
-    def matvec(self, /, x: Matrix[_Real]) -> Matrix[_SCT_co]: ...
+    def matvec(self, /, x: onp.Matrix[_Real]) -> onp.Matrix[_SCT_co]: ...
     @overload  # complex matrix
-    def matvec(self, /, x: Matrix[_Number]) -> Matrix[_SCT_co | np.complex128]: ...
+    def matvec(self, /, x: onp.Matrix[_Number]) -> onp.Matrix[_SCT_co | np.complex128]: ...
     @overload  # float array 2d
     def matvec(self, /, x: onp.ToFloatStrict2D) -> onp.Array2D[_SCT_co]: ...
     @overload  # complex array 1d
diff --git a/scipy-stubs/spatial/distance.pyi b/scipy-stubs/spatial/distance.pyi
index 42f17f31..2abc22c6 100644
--- a/scipy-stubs/spatial/distance.pyi
+++ b/scipy-stubs/spatial/distance.pyi
@@ -5,7 +5,6 @@ from typing_extensions import TypeVar, deprecated
 
 import numpy as np
 import optype.numpy as onp
-import optype.typing as opt
 from scipy._typing import ToRNG
 
 __all__ = [
@@ -193,11 +192,11 @@ def pdist(X: onp.ToComplex2D, metric: _MetricFunc, *, out: _ArrayT, **kwargs: ob
 
 #
 @overload  # 1-d int
-def squareform(X: Seq[opt.JustInt], force: _Force = "no", checks: bool = True) -> onp.Array2D[np.int_]: ...
+def squareform(X: onp.ToJustIntStrict1D, force: _Force = "no", checks: bool = True) -> onp.Array2D[np.int_]: ...
 @overload  # 1-d float
-def squareform(X: Seq[opt.Just[float]], force: _Force = "no", checks: bool = True) -> onp.Array2D[np.float64]: ...
+def squareform(X: onp.ToJustFloatStrict1D, force: _Force = "no", checks: bool = True) -> onp.Array2D[np.float64]: ...
 @overload  # 1-d complex
-def squareform(X: Seq[opt.Just[complex]], force: _Force = "no", checks: bool = True) -> onp.Array2D[np.complex128]: ...
+def squareform(X: onp.ToJustComplexStrict1D, force: _Force = "no", checks: bool = True) -> onp.Array2D[np.complex128]: ...
 @overload  # 1-d array-like
 def squareform(
     X: Seq[_NumberT] | onp.CanArray1D[_NumberT],
@@ -205,11 +204,11 @@ def squareform(
     checks: bool = True,
 ) -> onp.Array2D[_NumberT]: ...
 @overload  # 2-d int
-def squareform(X: _Seq2D[opt.JustInt], force: _Force = "no", checks: bool = True) -> onp.Array1D[np.int_]: ...
+def squareform(X: onp.ToJustIntStrict2D, force: _Force = "no", checks: bool = True) -> onp.Array1D[np.int_]: ...
 @overload  # 2-d float
-def squareform(X: _Seq2D[opt.Just[float]], force: _Force = "no", checks: bool = True) -> onp.Array1D[np.float64]: ...
+def squareform(X: onp.ToJustFloatStrict2D, force: _Force = "no", checks: bool = True) -> onp.Array1D[np.float64]: ...
 @overload  # 2-d complex
-def squareform(X: _Seq2D[opt.Just[complex]], force: _Force = "no", checks: bool = True) -> onp.Array1D[np.complex128]: ...
+def squareform(X: onp.ToJustComplexStrict2D, force: _Force = "no", checks: bool = True) -> onp.Array1D[np.complex128]: ...
 @overload  # 2-d array-like
 def squareform(
     X: _Seq2D[_NumberT] | Seq[onp.Array1D[_NumberT]] | onp.CanArray2D[_NumberT],
diff --git a/scipy-stubs/special/_multiufuncs.pyi b/scipy-stubs/special/_multiufuncs.pyi
index 69445b6c..a4bd6949 100644
--- a/scipy-stubs/special/_multiufuncs.pyi
+++ b/scipy-stubs/special/_multiufuncs.pyi
@@ -22,7 +22,7 @@ _RT = TypeVar("_RT")
 _UFuncT_co = TypeVar("_UFuncT_co", bound=Callable[..., object], default=Callable[..., Any], covariant=True)
 
 _Complex: TypeAlias = np.complex64 | np.complex128  # `clongdouble` isn't supported
-_ToJustComplex: TypeAlias = opt.Just[complex] | _Complex
+_ToJustComplex: TypeAlias = opt.JustComplex | _Complex
 _ToJustComplexND: TypeAlias = onp.CanArrayND[_Complex] | onp.SequenceND[onp.CanArrayND[_Complex]] | onp.SequenceND[_ToJustComplex]
 _ToJustComplex_D: TypeAlias = _ToJustComplex | _ToJustComplexND
 
diff --git a/scipy-stubs/special/_ufuncs.pyi b/scipy-stubs/special/_ufuncs.pyi
index cf81efcd..641b6f56 100644
--- a/scipy-stubs/special/_ufuncs.pyi
+++ b/scipy-stubs/special/_ufuncs.pyi
@@ -326,10 +326,10 @@ _CoFloat64ND: TypeAlias = onp.ArrayND[_CoFloat64]
 _CoComplex128ND: TypeAlias = onp.ArrayND[_CoComplex128]
 
 _SubFloat: TypeAlias = _Float16 | _CoInt  # anything "below" float32 | float64 that isn't float32 | float64
-_ToSubFloat: TypeAlias = opt.Just[float] | int | _SubFloat  # does not overlap with float32 | float64
-_ToSubFloatND: TypeAlias = _ToND[_SubFloat, opt.Just[float] | int]
+_ToSubFloat: TypeAlias = opt.JustFloat | int | _SubFloat  # does not overlap with float32 | float64
+_ToSubFloatND: TypeAlias = _ToND[_SubFloat, opt.JustFloat | int]
 
-_ToSubComplex: TypeAlias = opt.Just[complex] | _ToSubFloat  # does not overlap with complex64 | complex128
+_ToSubComplex: TypeAlias = opt.JustComplex | _ToSubFloat  # does not overlap with complex64 | complex128
 
 _CoT = TypeVar("_CoT", bound=np.generic)
 _ToT = TypeVar("_ToT")
@@ -854,11 +854,11 @@ class _UFunc11fc(_UFunc11[_NameT_co, _IdentityT_co], Generic[_NameT_co, _Identit
     def types(self, /) -> list[L["f->f", "d->d", "F->F", "D->D"]]: ...
     #
     @overload
-    def __call__(self, x: opt.Just[float] | opt.JustInt, /, out: _Out1 = None, **kw: Unpack[_KwBase]) -> _Float64: ...
+    def __call__(self, x: opt.JustFloat | opt.JustInt, /, out: _Out1 = None, **kw: Unpack[_KwBase]) -> _Float64: ...
     @overload
     def __call__(self, x: _ToSubFloat, /, out: _Out1 = None, **kw: Unpack[_Kw11fc]) -> _Float: ...
     @overload
-    def __call__(self, x: opt.Just[complex], /, out: _Out1 = None, **kw: Unpack[_KwBase]) -> _Complex128: ...
+    def __call__(self, x: opt.JustComplex, /, out: _Out1 = None, **kw: Unpack[_KwBase]) -> _Complex128: ...
     @overload
     def __call__(self, x: _ToSubComplex, /, out: _Out1 = None, **kw: Unpack[_Kw11fc]) -> _Inexact: ...
     @overload
diff --git a/scipy-stubs/stats/_distribution_infrastructure.pyi b/scipy-stubs/stats/_distribution_infrastructure.pyi
index c0994d63..212b3fe5 100644
--- a/scipy-stubs/stats/_distribution_infrastructure.pyi
+++ b/scipy-stubs/stats/_distribution_infrastructure.pyi
@@ -62,8 +62,7 @@ _ToJustIntMax2D: TypeAlias = onp.ToJustIntStrict2D | _ToJustIntMax1D
 _ToJustIntMax3D: TypeAlias = onp.ToJustIntStrict3D | _ToJustIntMax2D
 _ToJustIntMaxND: TypeAlias = onp.ToJustIntND | onp.ToJustInt
 
-_JustFloat: TypeAlias = opt.Just[float] | _Float
-_Null: TypeAlias = opt.Just[object]
+_Null: TypeAlias = opt.Just[object]  # type of `_null`
 _Axes: TypeAlias = object  # placeholder for `matplotlib.axes.Axes`
 
 _DomainRegion: TypeAlias = L["domain", "typical"]
@@ -75,7 +74,7 @@ _SMomentMethod: TypeAlias = L["formula", "general", "transform", "normalize", "c
 
 _ParamValues: TypeAlias = Mapping[str, _ToFloatMaxND]
 _ToDomain: TypeAlias = tuple[onp.ToFloat | str, onp.ToFloat | str]
-_ToTol: TypeAlias = opt.Just[float] | _Null
+_ToTol: TypeAlias = opt.JustFloat | _Null
 _DrawProportions: TypeAlias = tuple[onp.ToFloat, onp.ToFloat, onp.ToFloat, onp.ToFloat]
 _Elementwise: TypeAlias = Callable[[onp.ArrayND[np.float64]], onp.ArrayND[_FloatT]]
 
@@ -274,7 +273,9 @@ class ContinuousDistribution(_BaseDistribution[_FloatT_co, _ShapeT_co], Generic[
     def __sub__(self, lshift: onp.ToFloat, /) -> _LinDist[Self, _Float, _ShapeT_co]: ...
     @overload
     def __sub__(
-        self: _DistT0, lshift: onp.CanArrayND[_FloatT, _ShapeT], /
+        self: _DistT0,
+        lshift: onp.CanArrayND[_FloatT, _ShapeT],
+        /,
     ) -> _LinDist[_DistT0, _FloatT | _FloatT_co, _ShapeT]: ...
     @overload
     def __sub__(self: _DistT_1, lshift: onp.ToFloatStrict1D, /) -> _LinDist[_DistT_1, _Float, _1D]: ...
@@ -295,7 +296,9 @@ class ContinuousDistribution(_BaseDistribution[_FloatT_co, _ShapeT_co], Generic[
     def __mul__(self, scale: onp.ToFloat, /) -> _LinDist[Self, _Float, _ShapeT_co]: ...
     @overload
     def __mul__(
-        self: _DistT0, scale: onp.CanArrayND[_FloatT, _ShapeT], /
+        self: _DistT0,
+        scale: onp.CanArrayND[_FloatT, _ShapeT],
+        /,
     ) -> _LinDist[_DistT0, _FloatT | _FloatT_co, _ShapeT]: ...
     @overload
     def __mul__(self: _DistT_1, scale: onp.ToFloatStrict1D, /) -> _LinDist[_DistT_1, _Float, _1D]: ...
@@ -316,7 +319,9 @@ class ContinuousDistribution(_BaseDistribution[_FloatT_co, _ShapeT_co], Generic[
     def __truediv__(self, iscale: onp.ToFloat, /) -> _LinDist[Self, _Float, _ShapeT_co]: ...
     @overload
     def __truediv__(
-        self: _DistT0, iscale: onp.CanArrayND[_FloatT, _ShapeT], /
+        self: _DistT0,
+        iscale: onp.CanArrayND[_FloatT, _ShapeT],
+        /,
     ) -> _LinDist[_DistT0, _FloatT | _FloatT_co, _ShapeT]: ...
     @overload
     def __truediv__(self: _DistT_1, iscale: onp.ToFloatStrict1D, /) -> _LinDist[_DistT_1, _Float, _1D]: ...
@@ -339,7 +344,7 @@ class ContinuousDistribution(_BaseDistribution[_FloatT_co, _ShapeT_co], Generic[
         x: _PlotQuantity = "x",
         y: _PlotQuantity = "pdf",
         *,
-        t: tuple[_PlotQuantity, _JustFloat, _JustFloat] = ("cdf", 0.0005, 0.9995),
+        t: tuple[_PlotQuantity, onp.ToJustFloat, onp.ToJustFloat] = ("cdf", 0.0005, 0.9995),
         ax: _AxesT | None = None,
     ) -> _AxesT: ...
 
diff --git a/scipy-stubs/stats/_mstats_basic.pyi b/scipy-stubs/stats/_mstats_basic.pyi
index 80338603..120ab664 100644
--- a/scipy-stubs/stats/_mstats_basic.pyi
+++ b/scipy-stubs/stats/_mstats_basic.pyi
@@ -1,5 +1,5 @@
 from collections.abc import Callable
-from typing import Any, Final, Generic, Literal, NamedTuple, TypeAlias, TypedDict, overload, type_check_only
+from typing import Any, Concatenate, Final, Generic, Literal, NamedTuple, TypeAlias, TypedDict, overload, type_check_only
 from typing_extensions import Self, TypeVar
 
 import numpy as np
@@ -73,24 +73,27 @@ __all__ = [
     "winsorize",
 ]
 
-_ShapeT = TypeVar("_ShapeT", bound=tuple[int, ...])
-
 _SCT = TypeVar("_SCT", bound=np.generic, default=np.float64)
 _SCT_f = TypeVar("_SCT_f", bound=np.floating[Any], default=np.float64)
 _SCT_bifc = TypeVar("_SCT_bifc", bound=np.number[Any] | np.bool_, default=np.float64)
 _SCT_bifcmO = TypeVar("_SCT_bifcmO", bound=np.number[Any] | np.timedelta64 | np.bool_ | np.object_)
 
-_NDT_f_co = TypeVar("_NDT_f_co", covariant=True, bound=float | _MArrayND0[np.floating[Any]], default=_MArrayND[np.float64])
+_MArrayOrND: TypeAlias = _SCT | onp.MArray[_SCT]
+
+_NDT_f_co = TypeVar("_NDT_f_co", covariant=True, bound=float | _MArrayOrND[np.floating[Any]], default=onp.MArray[np.float64])
 _NDT_fc_co = TypeVar(
     "_NDT_fc_co",
     covariant=True,
-    bound=complex | _MArrayND0[np.inexact[Any]],
-    default=_MArrayND0[np.float64 | np.complex128],
+    bound=complex | _MArrayOrND[np.inexact[Any]],
+    default=_MArrayOrND[np.float64 | np.complex128],
 )
 
-_MArray: TypeAlias = np.ma.MaskedArray[_ShapeT, np.dtype[_SCT]]
-_MArrayND: TypeAlias = _MArray[tuple[int, ...], _SCT]
-_MArrayND0: TypeAlias = _SCT | _MArray[tuple[int, ...], _SCT]
+_KendallTauMethod: TypeAlias = Literal["auto", "asymptotic", "exact"]
+_TheilSlopesMethod: TypeAlias = Literal["joint", "separate"]
+_SiegelSlopesMethod: TypeAlias = Literal["hierarchical", "separate"]
+
+_KSMethod: TypeAlias = Literal["auto", "exact", "asymp"]
+_KTestMethod: TypeAlias = Literal[_KSMethod, "approx"]
 
 @type_check_only
 class _TestResult(NamedTuple, Generic[_NDT_f_co, _NDT_fc_co]):
@@ -100,14 +103,14 @@ class _TestResult(NamedTuple, Generic[_NDT_f_co, _NDT_fc_co]):
 _KendallTauSeasonalResult = TypedDict(
     "_KendallTauSeasonalResult",
     {
-        "seasonal tau": _MArrayND0[np.float64],
+        "seasonal tau": _MArrayOrND[np.float64],
         "global tau": np.float64,
         "global tau (alt)": np.float64,
         "seasonal p-value": onp.ArrayND[np.float64],
         "global p-value (indep)": np.float64,
         "global p-value (dep)": np.float64,
-        "chi2 total": _MArrayND[np.float64],
-        "chi2 trend": _MArrayND[np.float64],
+        "chi2 total": onp.MArray[np.float64],
+        "chi2 trend": onp.MArray[np.float64],
     },
 )
 
@@ -116,12 +119,12 @@ _KendallTauSeasonalResult = TypedDict(
 trimdoc: Final[str] = ...
 
 class ModeResult(NamedTuple):
-    mode: _MArrayND
-    count: _MArrayND  # type: ignore[assignment]  # pyright: ignore[reportIncompatibleMethodOverride]
+    mode: onp.MArray[np.float64]
+    count: onp.MArray[np.float64]  # type: ignore[assignment]  # pyright: ignore[reportIncompatibleMethodOverride]
 
 class DescribeResult(NamedTuple):
     nobs: np.int_ | onp.ArrayND[np.int_]
-    minmax: tuple[_MArrayND[np.floating[Any] | np.integer[Any]], _MArrayND[np.floating[Any] | np.integer[Any]]]
+    minmax: tuple[onp.MArray[np.floating[Any] | np.integer[Any]], onp.MArray[np.floating[Any] | np.integer[Any]]]
     mean: np.floating[Any]
     variance: np.floating[Any]
     skewness: np.floating[Any]
@@ -143,18 +146,18 @@ class KruskalResult(_TestResult[np.float64, np.float64]): ...
 class FriedmanchisquareResult(_TestResult[np.float64, np.float64]): ...
 class BrunnerMunzelResult(_TestResult[np.float64, np.float64]): ...
 
-class SenSeasonalSlopesResult(BaseBunch[_MArrayND[np.float64], np.float64]):
+class SenSeasonalSlopesResult(BaseBunch[onp.MArray[np.float64], np.float64]):
     def __new__(_cls, intra_slope: float, inter_slope: float) -> Self: ...
     def __init__(self, /, intra_slope: float, inter_slope: float) -> None: ...
     @property
-    def intra_slope(self, /) -> _MArrayND[np.float64]: ...
+    def intra_slope(self, /) -> onp.MArray[np.float64]: ...
     @property
     def inter_slope(self, /) -> float: ...
 
 # TODO(jorenham): Overloads for scalar vs. array
 # TODO(jorenham): Overloads for specific dtypes
 
-def argstoarray(*args: onp.ToFloatND) -> _MArrayND[np.float64]: ...
+def argstoarray(*args: onp.ToFloatND) -> onp.MArray[np.float64]: ...
 def find_repeats(arr: onp.ToFloatND) -> tuple[onp.ArrayND[np.float64], onp.ArrayND[np.intp]]: ...
 def count_tied_groups(x: onp.ToFloatND, use_missing: bool = False) -> dict[np.intp, np.intp]: ...
 def rankdata(data: onp.ToFloatND, axis: op.CanIndex | None = None, use_missing: bool = False) -> onp.ArrayND[np.float64]: ...
@@ -181,7 +184,7 @@ def kendalltau(
     y: onp.ToFloatND,
     use_ties: bool = True,
     use_missing: bool = False,
-    method: Literal["auto", "asymptotic", "exact"] = "auto",
+    method: _KendallTauMethod = "auto",
     alternative: Alternative = "two-sided",
 ) -> SignificanceResult: ...
 def kendalltau_seasonal(x: onp.ToFloatND) -> _KendallTauSeasonalResult: ...
@@ -190,13 +193,13 @@ def linregress(x: onp.ToFloatND, y: onp.ToFloatND | None = None) -> LinregressRe
 def theilslopes(
     y: onp.ToFloatND,
     x: onp.ToFloatND | None = None,
-    alpha: float = 0.95,
-    method: Literal["joint", "separate"] = "separate",
+    alpha: onp.ToJustFloat = 0.95,
+    method: _TheilSlopesMethod = "separate",
 ) -> TheilslopesResult: ...
 def siegelslopes(
     y: onp.ToFloatND,
     x: onp.ToFloatND | None = None,
-    method: Literal["hierarchical", "separate"] = "hierarchical",
+    method: _SiegelSlopesMethod = "hierarchical",
 ) -> SiegelslopesResult: ...
 def sen_seasonal_slopes(x: onp.ToFloatND) -> SenSeasonalSlopesResult: ...
 
@@ -224,25 +227,47 @@ def mannwhitneyu(x: onp.ToFloatND, y: onp.ToFloatND, use_continuity: op.CanBool
 def kruskal(arg0: onp.ToFloatND, arg1: onp.ToFloatND, /, *args: onp.ToFloatND) -> KruskalResult: ...
 
 #
+@overload
 def ks_1samp(
     x: onp.ToFloatND,
     cdf: str | Callable[[float], onp.ToFloat],
-    args: tuple[object, ...] = (),
+    args: tuple[()] = (),
+    alternative: Alternative = "two-sided",
+    method: _KSMethod = "auto",
+) -> KstestResult: ...
+@overload
+def ks_1samp(
+    x: onp.ToFloatND,
+    cdf: str | Callable[Concatenate[float, ...], onp.ToFloat],
+    args: tuple[object, ...],
     alternative: Alternative = "two-sided",
-    method: Literal["auto", "exact", "asymp"] = "auto",
+    method: _KSMethod = "auto",
 ) -> KstestResult: ...
+
+#
 def ks_2samp(
     data1: onp.ToFloatND,
     data2: onp.ToFloatND,
     alternative: Alternative = "two-sided",
-    method: Literal["auto", "exact", "asymp"] = "auto",
+    method: _KSMethod = "auto",
 ) -> KstestResult: ...
+
+#
+@overload
 def kstest(
     data1: onp.ToFloatND,
     data2: onp.ToFloatND | str | Callable[[float], onp.ToFloat],
-    args: tuple[object, ...] = (),
+    args: tuple[()] = (),
+    alternative: Alternative = "two-sided",
+    method: _KTestMethod = "auto",
+) -> KstestResult: ...
+@overload
+def kstest(
+    data1: onp.ToFloatND,
+    data2: Callable[Concatenate[float, ...], onp.ToFloat],
+    args: tuple[object, ...],
     alternative: Alternative = "two-sided",
-    method: Literal["auto", "exact", "approx", "asymp"] = "auto",
+    method: _KTestMethod = "auto",
 ) -> KstestResult: ...
 
 #
@@ -251,31 +276,31 @@ def trima(
     a: onp.SequenceND[bool],
     limits: tuple[onp.ToInt, onp.ToInt] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
-) -> _MArrayND[np.bool_]: ...
+) -> onp.MArray[np.bool_]: ...
 @overload
 def trima(
     a: onp.SequenceND[opt.JustInt],
     limits: tuple[onp.ToInt, onp.ToInt] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
-) -> _MArrayND[np.int_]: ...
+) -> onp.MArray[np.int_]: ...
 @overload
 def trima(
     a: onp.SequenceND[float],
     limits: tuple[onp.ToFloat, onp.ToFloat] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
-) -> _MArrayND[np.float64 | np.int_ | np.bool_]: ...
+) -> onp.MArray[np.float64 | np.int_ | np.bool_]: ...
 @overload
 def trima(
     a: onp.SequenceND[complex],
     limits: tuple[onp.ToComplex, onp.ToComplex] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
-) -> _MArrayND[np.complex128 | np.float64 | np.int_ | np.bool_]: ...
+) -> onp.MArray[np.complex128 | np.float64 | np.int_ | np.bool_]: ...
 @overload
 def trima(
     a: _ArrayLike[_SCT_bifc],
     limits: tuple[onp.ToComplex, onp.ToComplex] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
-) -> _MArrayND[_SCT_bifc]: ...
+) -> onp.MArray[_SCT_bifc]: ...
 
 #
 @overload
@@ -284,28 +309,28 @@ def trimr(
     limits: tuple[onp.ToFloat, onp.ToFloat] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.int_]: ...
+) -> onp.MArray[np.int_]: ...
 @overload
 def trimr(
     a: onp.SequenceND[float],
     limits: tuple[onp.ToFloat, onp.ToFloat] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.float64 | np.int_]: ...
+) -> onp.MArray[np.float64 | np.int_]: ...
 @overload
 def trimr(
     a: onp.SequenceND[complex],
     limits: tuple[onp.ToComplex, onp.ToComplex] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.complex128 | np.float64 | np.int_]: ...
+) -> onp.MArray[np.complex128 | np.float64 | np.int_]: ...
 @overload
 def trimr(
     a: _ArrayLike[_SCT_bifc],
     limits: tuple[onp.ToComplex, onp.ToComplex] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[_SCT_bifc]: ...
+) -> onp.MArray[_SCT_bifc]: ...
 
 #
 @overload
@@ -315,7 +340,7 @@ def trim(
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     relative: op.CanBool = False,
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.int_]: ...
+) -> onp.MArray[np.int_]: ...
 @overload
 def trim(
     a: onp.SequenceND[float],
@@ -323,7 +348,7 @@ def trim(
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     relative: op.CanBool = False,
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.float64 | np.int_]: ...
+) -> onp.MArray[np.float64 | np.int_]: ...
 @overload
 def trim(
     a: onp.SequenceND[complex],
@@ -331,7 +356,7 @@ def trim(
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     relative: op.CanBool = False,
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.complex128 | np.float64 | np.int_]: ...
+) -> onp.MArray[np.complex128 | np.float64 | np.int_]: ...
 @overload
 def trim(
     a: _ArrayLike[_SCT_bifc],
@@ -339,7 +364,7 @@ def trim(
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     relative: op.CanBool = False,
     axis: op.CanIndex | None = None,
-) -> _MArrayND[_SCT_bifc]: ...
+) -> onp.MArray[_SCT_bifc]: ...
 
 #
 @overload
@@ -348,28 +373,28 @@ def trimboth(
     proportiontocut: float | np.floating[Any] = 0.2,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.int_]: ...
+) -> onp.MArray[np.int_]: ...
 @overload
 def trimboth(
     data: onp.SequenceND[float],
     proportiontocut: float | np.floating[Any] = 0.2,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.float64 | np.int_]: ...
+) -> onp.MArray[np.float64 | np.int_]: ...
 @overload
 def trimboth(
     data: onp.SequenceND[complex],
     proportiontocut: float | np.floating[Any] = 0.2,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.complex128 | np.float64 | np.int_]: ...
+) -> onp.MArray[np.complex128 | np.float64 | np.int_]: ...
 @overload
 def trimboth(
     data: _ArrayLike[_SCT_bifc],
     proportiontocut: float | np.floating[Any] = 0.2,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[_SCT_bifc]: ...
+) -> onp.MArray[_SCT_bifc]: ...
 
 #
 @overload
@@ -379,7 +404,7 @@ def trimtail(
     tail: Literal["left", "right"] = "left",
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.int_]: ...
+) -> onp.MArray[np.int_]: ...
 @overload
 def trimtail(
     data: onp.SequenceND[float],
@@ -387,7 +412,7 @@ def trimtail(
     tail: Literal["left", "right"] = "left",
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.float64 | np.int_]: ...
+) -> onp.MArray[np.float64 | np.int_]: ...
 @overload
 def trimtail(
     data: onp.SequenceND[complex],
@@ -395,7 +420,7 @@ def trimtail(
     tail: Literal["left", "right"] = "left",
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[np.complex128 | np.float64 | np.int_]: ...
+) -> onp.MArray[np.complex128 | np.float64 | np.int_]: ...
 @overload
 def trimtail(
     data: _ArrayLike[_SCT_bifc],
@@ -403,7 +428,7 @@ def trimtail(
     tail: Literal["left", "right"] = "left",
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND[_SCT_bifc]: ...
+) -> onp.MArray[_SCT_bifc]: ...
 
 #
 @overload
@@ -413,7 +438,7 @@ def trimmed_mean(
     inclusive: tuple[op.CanBool, op.CanBool] = (1, 1),
     relative: op.CanBool = True,
     axis: op.CanIndex | None = None,
-) -> _MArrayND0[np.floating[Any]]: ...
+) -> _MArrayOrND[np.floating[Any]]: ...
 @overload
 def trimmed_mean(
     a: onp.ToComplexND,
@@ -421,7 +446,7 @@ def trimmed_mean(
     inclusive: tuple[op.CanBool, op.CanBool] = (1, 1),
     relative: op.CanBool = True,
     axis: op.CanIndex | None = None,
-) -> _MArrayND0[np.floating[Any] | np.complex128]: ...
+) -> _MArrayOrND[np.floating[Any] | np.complex128]: ...
 
 #
 def trimmed_var(
@@ -431,7 +456,7 @@ def trimmed_var(
     relative: op.CanBool = True,
     axis: op.CanIndex | None = None,
     ddof: onp.ToInt = 0,
-) -> _MArrayND0[np.float64]: ...
+) -> _MArrayOrND[np.float64]: ...
 
 #
 def trimmed_std(
@@ -441,7 +466,7 @@ def trimmed_std(
     relative: op.CanBool = True,
     axis: op.CanIndex | None = None,
     ddof: onp.ToInt = 0,
-) -> _MArrayND0[np.float64]: ...
+) -> _MArrayOrND[np.float64]: ...
 
 #
 def trimmed_stde(
@@ -449,7 +474,7 @@ def trimmed_stde(
     limits: tuple[onp.ToFloat, onp.ToFloat] = (0.1, 0.1),
     inclusive: tuple[op.CanBool, op.CanBool] = (1, 1),
     axis: op.CanIndex | None = None,
-) -> _MArrayND0[np.float64]: ...
+) -> _MArrayOrND[np.float64]: ...
 
 #
 @overload
@@ -458,23 +483,23 @@ def tmean(
     limits: tuple[onp.ToFloat, onp.ToFloat] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND0[np.floating[Any]]: ...
+) -> _MArrayOrND[np.floating[Any]]: ...
 @overload
 def tmean(
     a: onp.ToComplexND,
     limits: tuple[onp.ToComplex, onp.ToComplex] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = None,
-) -> _MArrayND0[np.inexact[Any]]: ...
+) -> _MArrayOrND[np.inexact[Any]]: ...
 
 #
 def tvar(
-    a: _MArrayND,
+    a: onp.MArray[np.floating[Any] | np.integer[Any]],
     limits: tuple[onp.ToFloat, onp.ToFloat] | None = None,
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = 0,
     ddof: onp.ToInt = 1,
-) -> _MArrayND0[np.floating[Any]]: ...
+) -> _MArrayOrND[np.floating[Any]]: ...
 
 #
 @overload
@@ -483,28 +508,28 @@ def tmin(
     lowerlimit: onp.ToFloat | None = None,
     axis: op.CanIndex | None = 0,
     inclusive: AnyBool = True,
-) -> _MArrayND0[np.int_]: ...
+) -> _MArrayOrND[np.int_]: ...
 @overload
 def tmin(
     a: onp.SequenceND[float],
     lowerlimit: onp.ToFloat | None = None,
     axis: op.CanIndex | None = 0,
     inclusive: AnyBool = True,
-) -> _MArrayND0[np.float64 | np.int_]: ...
+) -> _MArrayOrND[np.float64 | np.int_]: ...
 @overload
 def tmin(
     a: onp.SequenceND[complex],
     lowerlimit: onp.ToComplex | None = None,
     axis: op.CanIndex | None = 0,
     inclusive: AnyBool = True,
-) -> _MArrayND0[np.complex128 | np.float64 | np.int_]: ...
+) -> _MArrayOrND[np.complex128 | np.float64 | np.int_]: ...
 @overload
 def tmin(
     a: _ArrayLike[_SCT_bifc],
     lowerlimit: onp.ToComplex | None = None,
     axis: op.CanIndex | None = 0,
     inclusive: AnyBool = True,
-) -> _MArrayND0[_SCT_bifc]: ...
+) -> _MArrayOrND[_SCT_bifc]: ...
 
 #
 @overload
@@ -513,28 +538,28 @@ def tmax(
     upperlimit: onp.ToFloat | None = None,
     axis: op.CanIndex | None = 0,
     inclusive: AnyBool = True,
-) -> _MArrayND0[np.int_]: ...
+) -> _MArrayOrND[np.int_]: ...
 @overload
 def tmax(
     a: onp.SequenceND[float],
     upperlimit: onp.ToFloat | None = None,
     axis: op.CanIndex | None = 0,
     inclusive: AnyBool = True,
-) -> _MArrayND0[np.float64 | np.int_]: ...
+) -> _MArrayOrND[np.float64 | np.int_]: ...
 @overload
 def tmax(
     a: onp.SequenceND[complex],
     upperlimit: onp.ToComplex | None = None,
     axis: op.CanIndex | None = 0,
     inclusive: AnyBool = True,
-) -> _MArrayND0[np.complex128 | np.float64 | np.int_]: ...
+) -> _MArrayOrND[np.complex128 | np.float64 | np.int_]: ...
 @overload
 def tmax(
     a: _ArrayLike[_SCT_bifc],
     upperlimit: onp.ToComplex | None = None,
     axis: op.CanIndex | None = 0,
     inclusive: AnyBool = True,
-) -> _MArrayND0[_SCT_bifc]: ...
+) -> _MArrayOrND[_SCT_bifc]: ...
 
 #
 def tsem(
@@ -543,7 +568,7 @@ def tsem(
     inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
     axis: op.CanIndex | None = 0,
     ddof: onp.ToInt = 1,
-) -> _MArrayND0: ...
+) -> _MArrayOrND: ...
 
 #
 @overload
@@ -554,7 +579,7 @@ def winsorize(
     inplace: AnyBool = False,
     axis: op.CanIndex | None = None,
     nan_policy: NanPolicy = "propagate",
-) -> _MArrayND[np.int_]: ...
+) -> onp.MArray[np.int_]: ...
 @overload
 def winsorize(
     a: _ArrayLike[_SCT_f],
@@ -563,7 +588,7 @@ def winsorize(
     inplace: AnyBool = False,
     axis: op.CanIndex | None = None,
     nan_policy: NanPolicy = "propagate",
-) -> _MArrayND[_SCT_f]: ...
+) -> onp.MArray[_SCT_f]: ...
 @overload
 def winsorize(
     a: onp.ToFloatND,
@@ -572,7 +597,7 @@ def winsorize(
     inplace: AnyBool = False,
     axis: op.CanIndex | None = None,
     nan_policy: NanPolicy = "propagate",
-) -> _MArrayND[np.floating[Any] | np.int_]: ...
+) -> onp.MArray[np.floating[Any] | np.int_]: ...
 @overload
 def winsorize(
     a: onp.ToComplexND,
@@ -581,22 +606,22 @@ def winsorize(
     inplace: AnyBool = False,
     axis: op.CanIndex | None = None,
     nan_policy: NanPolicy = "propagate",
-) -> _MArrayND[np.complex128 | np.floating[Any] | np.int_]: ...
+) -> onp.MArray[np.complex128 | np.floating[Any] | np.int_]: ...
 
 # TODO(jorenham): Overloads for complex array-likes
 def moment(
     a: onp.ToFloatND,
     moment: onp.ToInt | onp.ToIntND = 1,
     axis: op.CanIndex | None = 0,
-) -> _MArrayND0[np.floating[Any]]: ...
-def variation(a: onp.ToFloatND, axis: op.CanIndex | None = 0, ddof: onp.ToInt = 0) -> _MArrayND0[np.floating[Any]]: ...
-def skew(a: onp.ToFloatND, axis: op.CanIndex | None = 0, bias: op.CanBool = True) -> _MArrayND0[np.floating[Any]]: ...
+) -> _MArrayOrND[np.floating[Any]]: ...
+def variation(a: onp.ToFloatND, axis: op.CanIndex | None = 0, ddof: onp.ToInt = 0) -> _MArrayOrND[np.floating[Any]]: ...
+def skew(a: onp.ToFloatND, axis: op.CanIndex | None = 0, bias: op.CanBool = True) -> _MArrayOrND[np.floating[Any]]: ...
 def kurtosis(
     a: onp.ToFloatND,
     axis: op.CanIndex | None = 0,
     fisher: op.CanBool = True,
     bias: op.CanBool = True,
-) -> _MArrayND0[np.floating[Any]]: ...
+) -> _MArrayOrND[np.floating[Any]]: ...
 def describe(
     a: onp.ToFloatND,
     axis: op.CanIndex | None = 0,
@@ -606,9 +631,9 @@ def describe(
 
 #
 @overload
-def stde_median(data: onp.ToFloatND, axis: op.CanIndex | None = None) -> _MArrayND0[np.floating[Any]]: ...
+def stde_median(data: onp.ToFloatND, axis: op.CanIndex | None = None) -> _MArrayOrND[np.floating[Any]]: ...
 @overload
-def stde_median(data: onp.ToComplexND, axis: op.CanIndex | None = None) -> _MArrayND0[np.inexact[Any]]: ...
+def stde_median(data: onp.ToComplexND, axis: op.CanIndex | None = None) -> _MArrayOrND[np.inexact[Any]]: ...
 
 #
 @overload
@@ -616,13 +641,13 @@ def skewtest(
     a: onp.ToFloatND,
     axis: op.CanIndex | None = 0,
     alternative: Alternative = "two-sided",
-) -> SkewtestResult[_MArrayND0[np.float64], _MArrayND0[np.float64]]: ...
+) -> SkewtestResult[_MArrayOrND[np.float64], _MArrayOrND[np.float64]]: ...
 @overload
 def skewtest(
     a: onp.ToComplexND,
     axis: op.CanIndex | None = 0,
     alternative: Alternative = "two-sided",
-) -> SkewtestResult[_MArrayND0[np.float64], _MArrayND0[np.float64 | np.complex128]]: ...
+) -> SkewtestResult[_MArrayOrND[np.float64], _MArrayOrND[np.float64 | np.complex128]]: ...
 
 #
 @overload
@@ -630,16 +655,16 @@ def kurtosistest(
     a: onp.ToFloatND,
     axis: op.CanIndex | None = 0,
     alternative: Alternative = "two-sided",
-) -> KurtosistestResult[_MArrayND0[np.float64], _MArrayND0[np.float64]]: ...
+) -> KurtosistestResult[_MArrayOrND[np.float64], _MArrayOrND[np.float64]]: ...
 @overload
 def kurtosistest(
     a: onp.ToComplexND,
     axis: op.CanIndex | None = 0,
     alternative: Alternative = "two-sided",
-) -> KurtosistestResult[_MArrayND0[np.float64], _MArrayND0[np.float64 | np.complex128]]: ...
+) -> KurtosistestResult[_MArrayOrND[np.float64], _MArrayOrND[np.float64 | np.complex128]]: ...
 
 #
-def normaltest(a: onp.ToFloatND, axis: op.CanIndex | None = 0) -> NormaltestResult[_MArrayND0[np.float64]]: ...
+def normaltest(a: onp.ToFloatND, axis: op.CanIndex | None = 0) -> NormaltestResult[_MArrayOrND[np.float64]]: ...
 
 #
 def mquantiles(
@@ -649,7 +674,7 @@ def mquantiles(
     betap: onp.ToFloat = 0.4,
     axis: op.CanIndex | None = None,
     limit: tuple[onp.ToFloat, onp.ToFloat] | tuple[()] = (),
-) -> _MArrayND: ...
+) -> onp.MArray[np.float64]: ...
 
 #
 def scoreatpercentile(
@@ -658,16 +683,16 @@ def scoreatpercentile(
     limit: tuple[onp.ToFloat, onp.ToFloat] | tuple[()] = (),
     alphap: onp.ToFloat = 0.4,
     betap: onp.ToFloat = 0.4,
-) -> _MArrayND: ...
+) -> onp.MArray[np.float64]: ...
 
 #
-def plotting_positions(data: onp.ToFloatND, alpha: onp.ToFloat = 0.4, beta: onp.ToFloat = 0.4) -> _MArrayND: ...
+def plotting_positions(data: onp.ToFloatND, alpha: onp.ToFloat = 0.4, beta: onp.ToFloat = 0.4) -> onp.MArray[np.float64]: ...
 
 #
-def obrientransform(arg0: onp.ToFloatND, /, *args: onp.ToFloatND) -> _MArrayND: ...
+def obrientransform(arg0: onp.ToFloatND, /, *args: onp.ToFloatND) -> onp.MArray[np.float64]: ...
 
 #
-def sem(a: onp.ToFloatND, axis: op.CanIndex | None = 0, ddof: onp.ToInt = 1) -> np.float64 | _MArrayND: ...
+def sem(a: onp.ToFloatND, axis: op.CanIndex | None = 0, ddof: onp.ToInt = 1) -> np.float64 | onp.MArray[np.float64]: ...
 
 #
 def f_oneway(arg0: onp.ToFloatND, arg1: onp.ToFloatND, /, *args: onp.ToFloatND) -> F_onewayResult: ...
diff --git a/scipy-stubs/stats/_mstats_extras.pyi b/scipy-stubs/stats/_mstats_extras.pyi
index 01e86557..c3afbe84 100644
--- a/scipy-stubs/stats/_mstats_extras.pyi
+++ b/scipy-stubs/stats/_mstats_extras.pyi
@@ -1,4 +1,5 @@
-from typing import Any, Literal, overload
+from typing import Literal, TypeAlias, overload
+from typing_extensions import TypeVar
 
 import numpy as np
 import optype as op
@@ -17,113 +18,97 @@ __all__ = [
     "trimmed_mean_ci",
 ]
 
+_T = TypeVar("_T")
+_Tuple2: TypeAlias = tuple[_T, _T]
+_Falsy: TypeAlias = Literal[False, 0]
+_Truthy: TypeAlias = Literal[True, 1]
+_FloatND: TypeAlias = onp.ArrayND[np.float64]
+
+_ToProb: TypeAlias = onp.ToFloat | onp.ToFloatND
+_ToAxis: TypeAlias = op.CanIndex | None
+
+###
+
 @overload
 def hdquantiles(
-    data: onp.ToFloatND,
-    prob: onp.ToFloat | onp.ToFloatND = [0.25, 0.5, 0.75],
-    axis: op.CanIndex | None = None,
-    var: Literal[0, False] = False,
-) -> np.ma.MaskedArray[onp.AtLeast1D, np.dtype[np.float64]]: ...
+    data: onp.ToFloat1D,
+    prob: _ToProb = [0.25, 0.5, 0.75],
+    axis: _ToAxis = None,
+    var: _Falsy = False,
+) -> onp.MArray1D[np.float64]: ...
 @overload
 def hdquantiles(
-    data: onp.ToFloatND,
-    prob: onp.ToFloat | onp.ToFloatND,
-    axis: op.CanIndex | None,
-    var: Literal[1, True],
-) -> np.ma.MaskedArray[onp.AtLeast2D, np.dtype[np.float64]]: ...
+    data: onp.ToFloat1D,
+    prob: _ToProb,
+    axis: _ToAxis,
+    var: _Truthy,
+) -> onp.MArray2D[np.float64]: ...
 @overload
 def hdquantiles(
-    data: onp.ToFloatND,
-    prob: onp.ToFloat | onp.ToFloatND = [0.25, 0.5, 0.75],
-    axis: op.CanIndex | None = None,
+    data: onp.ToFloat1D,
+    prob: _ToProb = [0.25, 0.5, 0.75],
+    axis: _ToAxis = None,
     *,
-    var: Literal[1, True],
-) -> np.ma.MaskedArray[onp.AtLeast2D, np.dtype[np.float64]]: ...
+    var: _Truthy,
+) -> onp.MArray2D[np.float64]: ...
+@overload
+def hdquantiles(
+    data: onp.ToFloatND,
+    prob: _ToProb = [0.25, 0.5, 0.75],
+    axis: _ToAxis = None,
+    var: bool = False,
+) -> onp.MArray[np.float64]: ...
 
 #
 @overload
-def hdmedian(
-    data: onp.ToFloatND,
-    axis: op.CanIndex | None = -1,
-    var: Literal[0, False] = False,
-) -> np.ma.MaskedArray[onp.AtLeast0D, np.dtype[np.float64]]: ...
+def hdmedian(data: onp.ToFloatND, axis: _ToAxis = -1, var: _Falsy = False) -> onp.MArray[np.float64]: ...
 @overload
-def hdmedian(
-    data: onp.ToFloatND,
-    axis: op.CanIndex | None,
-    var: Literal[1, True],
-) -> np.ma.MaskedArray[onp.AtLeast1D, np.dtype[np.float64]]: ...
+def hdmedian(data: onp.ToFloatND, axis: _ToAxis, var: _Truthy) -> onp.MArray[np.float64]: ...
 @overload
-def hdmedian(
-    data: onp.ToFloatND,
-    axis: op.CanIndex | None = -1,
-    *,
-    var: Literal[1, True],
-) -> np.ma.MaskedArray[onp.AtLeast1D, np.dtype[np.float64]]: ...
+def hdmedian(data: onp.ToFloatND, axis: _ToAxis = -1, *, var: _Truthy) -> onp.MArray[np.float64]: ...
 
 #
-def hdquantiles_sd(
-    data: onp.ToFloatND,
-    prob: onp.ToFloat | onp.ToFloatND = (0.25, 0.5, 0.75),
-    axis: op.CanIndex | None = None,
-) -> np.ma.MaskedArray[onp.AtLeast1D, np.dtype[np.float64]]: ...
+def hdquantiles_sd(data: onp.ToFloatND, prob: _ToProb = (0.25, 0.5, 0.75), axis: _ToAxis = None) -> onp.MArray[np.float64]: ...
 
 #
 def trimmed_mean_ci(
     data: onp.ToFloatND,
-    limits: tuple[onp.ToFloat, onp.ToFloat] | None = (0.2, 0.2),
-    inclusive: tuple[op.CanBool, op.CanBool] = (True, True),
-    alpha: float | np.floating[Any] = 0.05,
-    axis: op.CanIndex | None = None,
-) -> onp.ArrayND[np.float64]: ...
+    limits: _Tuple2[onp.ToFloat] | None = (0.2, 0.2),
+    inclusive: _Tuple2[op.CanBool] = (True, True),
+    alpha: onp.ToJustFloat = 0.05,
+    axis: _ToAxis = None,
+) -> _FloatND: ...
 
 #
-def mjci(
-    data: onp.ToFloatND,
-    prob: onp.ToFloat | onp.ToFloatND = (0.25, 0.5, 0.75),
-    axis: op.CanIndex | None = None,
-) -> onp.ArrayND[np.float64]: ...
+def mjci(data: onp.ToFloatND, prob: _ToProb = (0.25, 0.5, 0.75), axis: _ToAxis = None) -> _FloatND: ...
 
 #
 def mquantiles_cimj(
     data: onp.ToFloatND,
-    prob: onp.ToFloat | onp.ToFloatND = (0.25, 0.5, 0.75),
-    alpha: float | np.floating[Any] = 0.05,
-    axis: op.CanIndex | None = None,
-) -> tuple[onp.ArrayND[np.float64], onp.ArrayND[np.float64]]: ...
+    prob: _ToProb = (0.25, 0.5, 0.75),
+    alpha: onp.ToJustFloat = 0.05,
+    axis: _ToAxis = None,
+) -> _Tuple2[_FloatND]: ...
 
 #
 @overload
-def median_cihs(
-    data: onp.ToFloatND,
-    alpha: float | np.floating[Any] = 0.05,
-    axis: None = None,
-) -> tuple[np.float64, np.float64]: ...
+def median_cihs(data: onp.ToFloatND, alpha: onp.ToJustFloat = 0.05, axis: None = None) -> _Tuple2[np.float64]: ...
 @overload
-def median_cihs(
-    data: onp.ToFloatND,
-    alpha: float | np.floating[Any],
-    axis: op.CanIndex,
-) -> tuple[np.float64 | onp.ArrayND[np.float64], np.float64 | onp.ArrayND[np.float64]]: ...
+def median_cihs(data: onp.ToFloatND, alpha: onp.ToJustFloat, axis: op.CanIndex) -> _Tuple2[np.float64 | _FloatND]: ...
 @overload
-def median_cihs(
-    data: onp.ToFloatND,
-    alpha: float | np.floating[Any] = 0.05,
-    *,
-    axis: op.CanIndex,
-) -> tuple[np.float64 | onp.ArrayND[np.float64], np.float64 | onp.ArrayND[np.float64]]: ...
+def median_cihs(data: onp.ToFloatND, alpha: onp.ToJustFloat = 0.05, *, axis: op.CanIndex) -> _Tuple2[np.float64 | _FloatND]: ...
 
 #
 @overload
 def compare_medians_ms(group_1: onp.ToFloatND, group_2: onp.ToFloatND, axis: None = None) -> np.float64: ...
 @overload
-def compare_medians_ms(group_1: onp.ToFloatND, group_2: onp.ToFloatND, axis: op.CanIndex) -> onp.ArrayND[np.float64]: ...
+def compare_medians_ms(group_1: onp.ToFloatND, group_2: onp.ToFloatND, axis: op.CanIndex) -> _FloatND: ...
 
 #
 @overload
 def idealfourths(data: onp.ToFloatND, axis: None = None) -> list[np.float64]: ...
 @overload
-def idealfourths(data: onp.ToFloatND, axis: op.CanIndex) -> np.ma.MaskedArray[onp.AtLeast1D, np.dtype[np.float64]]: ...
+def idealfourths(data: onp.ToFloatND, axis: op.CanIndex) -> onp.MArray[np.float64]: ...
 
 #
 def rsh(data: onp.ToFloatND, points: onp.ToFloatND | None = None) -> np.float64: ...
diff --git a/scipy-stubs/stats/_qmc.pyi b/scipy-stubs/stats/_qmc.pyi
index 0d682678..a4d2cc1a 100644
--- a/scipy-stubs/stats/_qmc.pyi
+++ b/scipy-stubs/stats/_qmc.pyi
@@ -1,6 +1,6 @@
 import abc
 import numbers
-from collections.abc import Callable, Mapping, Sequence
+from collections.abc import Callable, Mapping
 from typing import Any, ClassVar, Concatenate, Final, Generic, Literal, Protocol, TypeAlias, overload, type_check_only
 from typing_extensions import Self, TypeVar, override
 
@@ -34,10 +34,6 @@ _EngineT_co = TypeVar("_EngineT_co", bound=QMCEngine[np.inexact[Any]], default=S
 
 _Real: TypeAlias = np.floating[Any] | np.integer[Any]
 
-_ToJustFloat: TypeAlias = opt.Just[float] | np.floating[Any]
-_ToJustFloat1D: TypeAlias = onp.CanArrayND[np.floating[Any]] | Sequence[_ToJustFloat]
-_ToJustFloat2D: TypeAlias = onp.CanArrayND[np.floating[Any]] | Sequence[_ToJustFloat1D]
-
 _MethodQMC: TypeAlias = Literal["random-cd", "lloyd"]
 _MethodDisc: TypeAlias = Literal["CD", "WD", "MD", "L2-star"]
 _MethodDist: TypeAlias = Literal["mindist", "mst"]
@@ -257,7 +253,7 @@ class MultinomialQMC(_QMCDistribution[_EngineT_co], Generic[_EngineT_co]):
     def __init__(
         self: MultinomialQMC[Sobol],
         /,
-        pvals: _ToJustFloat | _ToJustFloat1D,
+        pvals: onp.ToJustFloat | onp.ToJustFloat1D,
         n_trials: onp.ToJustInt,
         *,
         engine: None = None,
@@ -267,7 +263,7 @@ class MultinomialQMC(_QMCDistribution[_EngineT_co], Generic[_EngineT_co]):
     def __init__(
         self,
         /,
-        pvals: _ToJustFloat | _ToJustFloat1D,
+        pvals: onp.ToJustFloat | onp.ToJustFloat1D,
         n_trials: onp.ToJustInt,
         *,
         engine: _EngineT_co,
@@ -300,29 +296,29 @@ def discrepancy(
 
 #
 def geometric_discrepancy(
-    sample: _ToJustFloat2D,
+    sample: onp.ToJustFloat2D,
     method: _MethodDist = "mindist",
     metric: _Metric = "euclidean",
 ) -> float | np.float64: ...
 
 #
-def update_discrepancy(x_new: _ToJustFloat1D, sample: _ToJustFloat2D, initial_disc: opt.AnyFloat) -> float: ...
+def update_discrepancy(x_new: onp.ToJustFloat1D, sample: onp.ToJustFloat2D, initial_disc: opt.AnyFloat) -> float: ...
 def primes_from_2_to(n: onp.ToInt) -> onp.Array1D[np.int_]: ...
 def n_primes(n: onp.ToInt) -> list[int] | onp.Array1D[np.int_]: ...
 
 #
 def _select_optimizer(optimization: _MethodQMC | None, config: Mapping[str, object]) -> _Optimizer | None: ...
 def _random_cd(best_sample: _FloatArrayT, n_iters: onp.ToInt, n_nochange: onp.ToInt, rng: RNG) -> _FloatArrayT: ...
-def _l1_norm(sample: _ToJustFloat2D) -> float | np.float64: ...
+def _l1_norm(sample: onp.ToJustFloat2D) -> float | np.float64: ...
 def _lloyd_iteration(sample: _FloatArrayT, decay: onp.ToFloat, qhull_options: str | None) -> _FloatArrayT: ...
 def _lloyd_centroidal_voronoi_tessellation(
-    sample: _ToJustFloat2D,
+    sample: onp.ToJustFloat2D,
     *,
     tol: onp.ToFloat = 1e-5,
     maxiter: onp.ToJustInt = 10,
     qhull_options: str | None = None,
 ) -> onp.Array2D[np.float64]: ...
-def _ensure_in_unit_hypercube(sample: _ToJustFloat2D) -> onp.Array2D[np.float64]: ...
+def _ensure_in_unit_hypercube(sample: onp.ToJustFloat2D) -> onp.Array2D[np.float64]: ...
 
 #
 @overload
diff --git a/tests/special/test_ufuncs.pyi b/tests/special/test_ufuncs.pyi
index ac09f3ac..f3a0745f 100644
--- a/tests/special/test_ufuncs.pyi
+++ b/tests/special/test_ufuncs.pyi
@@ -52,7 +52,7 @@ assert_type(sp.cbrt(_f4), np.float32)
 assert_type(sp.cbrt(_f4_nd), _Float32ND)
 assert_type(sp.cbrt(_f8), np.float64)
 assert_type(sp.cbrt(_f8_nd), _Float64ND)
-sp.cbrt(_c16)  # type:ignore[arg-type]  # pyright: ignore[reportArgumentType, reportCallIssue]
+sp.cbrt(_c16)  # type:ignore[call-overload]  # pyright: ignore[reportArgumentType, reportCallIssue]
 sp.cbrt(_c16_nd)  # type:ignore[arg-type]  # pyright: ignore[reportArgumentType, reportCallIssue]
 assert_type(sp.cbrt(False), np.float64)
 assert_type(sp.cbrt([False]), _Float64ND)
@@ -60,7 +60,7 @@ assert_type(sp.cbrt(0), np.float64)
 assert_type(sp.cbrt([0]), _Float64ND)
 assert_type(sp.cbrt(0.0), np.float64)
 assert_type(sp.cbrt([0.0]), _Float64ND)
-sp.cbrt(0j)  # type:ignore[arg-type]  # pyright: ignore[reportArgumentType, reportCallIssue]
+sp.cbrt(0j)  # type:ignore[call-overload]  # pyright: ignore[reportArgumentType, reportCallIssue]
 sp.cbrt([0j])  # type:ignore[arg-type]  # pyright: ignore[reportArgumentType, reportCallIssue]
 assert_type(sp.cbrt.at(_b1_nd, _i), None)
 assert_type(sp.cbrt.at(_f8_nd, _i), None)
@@ -74,7 +74,7 @@ assert_type(sp.logit(_f4), np.float32)
 assert_type(sp.logit(_f4_nd), _Float32ND)
 assert_type(sp.logit(_f8), np.float64)
 assert_type(sp.logit(_f8_nd), _Float64ND)
-sp.logit(_c16)  # type:ignore[arg-type]  # pyright: ignore[reportArgumentType, reportCallIssue]
+sp.logit(_c16)  # type:ignore[call-overload]  # pyright: ignore[reportArgumentType, reportCallIssue]
 sp.logit(_c16_nd)  # type:ignore[arg-type]  # pyright: ignore[reportArgumentType, reportCallIssue]
 assert_type(sp.logit(0), np.float64)
 assert_type(sp.logit([0]), _Float64ND)