Repository: scipy/scipy
File: scipy/stats/_morestats.py (Python)

Code before:
    >>> rng = np.random.default_rng()
    >>> x = stats.norm.rvs(loc=5, scale=3, size=100, random_state=rng)
    >>> shapiro_test = stats.shapiro(x)
    >>> shapiro_test
    ShapiroResult(statistic=0.9813305735588074, pvalue=0.16855233907699585)
    >>> shapiro_test.statistic
    0.9813305735588074
    >>> shapiro_test.pvalue
    0.16855233907699585

    For a more detailed example, see :ref:`hypothesis_shapiro`.
    """
    # `x` is a an array and axis=-1 due to _axis_nan_policy decorator
    xp = array_namespace(x)

    N = x.shape[-1]
    if N < 3:
        raise ValueError("Data must be at least length 3.")

    y = xp.sort(x, axis=-1)
    y -= x[..., N//2:N//2+1]  # subtract the median (or a nearby value); see gh-15777

    y = xp_promote(y, force_floating=True, xp=xp)
    w, pw = _swilk(y, xp=xp)
    if N > 5000:
        warnings.warn("scipy.stats.shapiro: For N > 5000, computed p-value "
                      f"may not be accurate. Current N is {N}.",
                      stacklevel=2)

    # `w` and `pw` are always Python floats, which are double precision.
    # We want to ensure that they are NumPy floats, so until dtypes are
    # respected, we can explicitly convert each to float64 (faster than
    # `np.array([w, pw])`).
    return ShapiroResult(w[()], pw[()])


def _swilk_w(y, a, *, xp):
    # calculate Shapiro-Wilk statistic given sorted sample and weights
    # Follows [4] Section 2.1
    num = xp.vecdot(a, y, axis=-1) ** 2
    y_ = _demean(y, mean=xp.mean(y, axis=-1, keepdims=True), axis=-1, xp=xp)
    den = xp.vecdot(y_, y_, axis=-1)
    return num / den


def _swilk(y, *, xp):
    # calculate Shapiro-Wilk statistic and p-value from sorted sample
    n = y.shape[-1]

    if n == 3:
        # [2] Table 5 gives the first four digits

Review comment (suggestion):
`N//2:N//2+1` preserves the dimensionality / broadcastability.

Code after:
    Examples
    --------

    >>> import numpy as np
    >>> from scipy import stats
    >>> rng = np.random.default_rng()
    >>> x = stats.norm.rvs(loc=5, scale=3, size=100, random_state=rng)
    >>> shapiro_test = stats.shapiro(x)
    >>> shapiro_test
    ShapiroResult(statistic=0.9813305735588074, pvalue=0.16855233907699585)
    >>> shapiro_test.statistic
    0.9813305735588074
    >>> shapiro_test.pvalue
    0.16855233907699585

    For a more detailed example, see :ref:`hypothesis_shapiro`.
    """
    # `x` is an array and axis=-1 due to _axis_nan_policy decorator
    xp = array_namespace(x)

    N = x.shape[-1]
    if N < 3:
        raise ValueError("Data must be at least length 3.")

    y = xp.sort(x, axis=-1)
    y -= x[..., N//2:N//2+1]  # subtract the median (or a nearby value); see gh-15777

    y = xp_promote(y, force_floating=True, xp=xp)
    w, pw = _swilk(y, xp=xp)
    if N > 5000:
        warnings.warn("scipy.stats.shapiro: For N > 5000, computed p-value "
                      f"may not be accurate. Current N is {N}.",
                      stacklevel=2)

    # `w` and `pw` are always Python floats, which are double precision.
    # We want to ensure that they are NumPy floats, so until dtypes are
    # respected, we can explicitly convert each to float64 (faster than
    # `np.array([w, pw])`).
    return ShapiroResult(w[()], pw[()])


def _swilk_w(y, a, *, xp):
    # calculate Shapiro-Wilk statistic given sorted sample and weights
    # Follows [4] Section 2.1
    num = xp.vecdot(a, y, axis=-1) ** 2
    y_ = _demean(y, mean=xp.mean(y, axis=-1, keepdims=True), axis=-1, xp=xp)
    den = xp.vecdot(y_, y_, axis=-1)
    return num / den


def _swilk(y, *, xp):