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

Code before:
    # 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(np.float64(w), np.float64(pw))


def _swilk_w(y, a):
    # calculate Shapiro-Wilk statistic given sorted sample and weights
    # Follows [4] Section 2.1
    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)


Review comment (suggestion):
Until we are using `array_namespace`...
```suggestion
        pvalue = 1. - 6/np.pi * np.arccos(np.sqrt(W))
```

Code after:


@xp_capabilities()
@_axis_nan_policy_factory(ShapiroResult, n_samples=1, too_small=2, default_axis=None)
def shapiro(x, *, axis=None):
    r"""Perform the Shapiro-Wilk test for normality.

    The Shapiro-Wilk test tests the null hypothesis that the
    data was drawn from a normal distribution.

    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