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

Code before:
    num = np.vecdot(a, y, axis=-1) ** 2
    y_ = _demean(y, mean=np.mean(y, axis=-1, keepdims=True), axis=-1, xp=np)
    den = np.vecdot(y_, y_, axis=-1)
    return num / den


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

    # Follows [4] section 2.2
    # could calculate half the coefficients and get the rest by antisymmetry
    i = np.arange(1, n + 1)
    m = special.ndtri((i - 3 / 8) / (n + 1 / 4))
    u = n**(-0.5)
    mTm = np.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)
    anm1 = (cnm1 + 0.042981*u - 0.293762*u**2
            - 1.752461*u**3 + 5.682633*u**4 - 3.582633*u**5)
    phi = ((mTm - 2*mn**2) / (1 - 2*an**2) if n <= 5
           else (mTm - 2*mn**2 - 2*mnm1**2) / (1 - 2*an**2 - 2*anm1**2))
    a = phi**(-0.5) * m
    if n > 5:
        a[..., -2] = anm1
        a[..., 1] = -anm1
    a[..., -1] = an
    a[..., 0] = -an
    W = _swilk_w(y, a)

    # Follows [4] Table 1
    if n <= 11:
        u = n
        gamma = -2.273 + 0.459*u
        mu = 0.5440 - 0.39978*u + 0.025054*u**2 - 0.0006714*u**3
        sigma = np.exp(1.3822 - 0.77857*u + 0.062767*u**2 - 0.0020322*u**3)
        gW = -np.log(gamma - np.log(1 - W))
    else:
        u = np.log(n)

Review comment (suggestion):
<img width="872" height="128" alt="image" src="https://github.com/user-attachments/assets/9db2cdc1-8c63-482f-8e49-6e93c0ab6025" />


Code after:

    Parameters
    ----------
    x : array_like
        Array of sample data. Must contain at least three observations.
    axis : int or tuple of ints, default: 0
        If an int or tuple of ints, the axis or axes of the input along which
        to compute the statistic. The statistic of each axis-slice (e.g. row)
        of the input will appear in a corresponding element of the output.
        If ``None``, the input will be raveled before computing the statistic.

    Returns
    -------
    statistic : float
        The test statistic.
    p-value : float
        The p-value for the hypothesis test.

    See Also
    --------
    anderson : The Anderson-Darling test for normality
    kstest : The Kolmogorov-Smirnov test for goodness of fit.
    :ref:`hypothesis_shapiro` : Extended example

    Notes
    -----
    The algorithm used is described in [4]_, but censoring parameters as
    described are not implemented. For N > 5000 the W test statistic is
    accurate, but the p-value may not be.

    References
    ----------
    .. [1] https://www.itl.nist.gov/div898/handbook/prc/section2/prc213.htm
           :doi:`10.18434/M32189`
    .. [2] Shapiro, S. S. & Wilk, M.B, "An analysis of variance test for
           normality (complete samples)", Biometrika, 1965, Vol. 52,
           pp. 591-611, :doi:`10.2307/2333709`
    .. [3] Razali, N. M. & Wah, Y. B., "Power comparisons of Shapiro-Wilk,
           Kolmogorov-Smirnov, Lilliefors and Anderson-Darling tests", Journal
           of Statistical Modeling and Analytics, 2011, Vol. 2, pp. 21-33.
    .. [4] Royston, P. "A toolkit for testing for non-normality in complete and
           censored samples." Journal of the Royal Statistical Society: Series D
           (The Statistician) 42.1 (1993): 37-43.

    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)