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

Code before:
    if N < 3:
        raise ValueError("Data must be at least length 3.")

    y = sort(x)
    y -= x[N//2]  # subtract the median (or a nearby value); see gh-15777

    w, pw = _swilk(y)
    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(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

Review comment (suggestion):
<img width="864" height="164" alt="image" src="https://github.com/user-attachments/assets/141e6296-3e4c-4998-b56c-96103d906b93" />


Code after:
    Determine and plot the optimal ``lmbda`` to transform ``x`` and plot it in
    the same plot:

    >>> _, maxlog = stats.yeojohnson(x)
    >>> ax.axvline(maxlog, color='r')

    >>> plt.show()

    """
    return _normplot('yeojohnson', x, la, lb, plot, N)


ShapiroResult = namedtuple('ShapiroResult', ('statistic', 'pvalue'))


@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.