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

Code before:

    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
        c = math.sqrt(2) / 2
        a = xp.asarray([-c, 0, c])
        # [2] Corollary 4; discussed in https://github.com/scipy/scipy/issues/18322
        W = xp.clip(_swilk_w(y, a, xp=xp), 0.75, 1.)
        pvalue = xp.clip(1. - 6/np.pi * xp.acos(xp.sqrt(W)), 0., 1.)
        return W, pvalue

    # Follows [4] section 2.2
    # could calculate half the coefficients and get the rest by antisymmetry
    i = xp.arange(1, n + 1, dtype=y.dtype)
    m = special.ndtri((i - 3 / 8) / (n + 1 / 4))
    u = n**(-0.5)
    mTm = xp.vecdot(m, m)
    c = mTm**(-0.5) * m
    mnm1, mn = m[..., -2], m[..., -1]
    cnm1, cn = c[..., -2], c[..., -1]
    an = (cn + 0.221157*u - 0.147981*u**2
          - 2.071190*u**3 + 4.434685*u**4 - 2.706056*u**5)

Review comment (suggestion):
<img width="1134" height="175" alt="image" src="https://github.com/user-attachments/assets/91cf9f64-bdb0-4496-8b78-c83100f4f64e" />


Code after:
    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
        c = math.sqrt(2) / 2
        a = xp.asarray([-c, 0, c])
        # [2] Corollary 4; discussed in https://github.com/scipy/scipy/issues/18322
        W = xp.clip(_swilk_w(y, a, xp=xp), 0.75, 1.)
        pvalue = xp.clip(1. - 6/np.pi * xp.acos(xp.sqrt(W)), 0., 1.)
        return W, pvalue

    # Follows [4] section 2.2
    # could calculate half the coefficients and get the rest by antisymmetry
    i = xp.arange(1, n + 1, dtype=y.dtype)
    m = special.ndtri((i - 3 / 8) / (n + 1 / 4))
    u = n**(-0.5)
    mTm = xp.vecdot(m, m)