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

Code before:
    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)
        mu = -1.5861 - 0.31082*u - 0.083751*u**2 + 0.0038915*u**3
        sigma = np.exp(-0.4803 - 0.082676*u + 0.0030302*u**2)
        gW = np.log(1 - W)

    z = (gW - mu) / sigma
    pvalue = special.ndtr(-z)
    return W, pvalue


# Values from [8]
_Avals_norm = array([0.561, 0.631, 0.752, 0.873, 1.035])
_Avals_expon = array([0.916, 1.062, 1.321, 1.591, 1.959])
# From Stephens, M A, "Goodness of Fit for the Extreme Value Distribution",
#             Biometrika, Vol. 64, Issue 3, Dec. 1977, pp 583-588.
_Avals_gumbel = array([0.474, 0.637, 0.757, 0.877, 1.038])
# From Stephens, M A, "Tests of Fit for the Logistic Distribution Based
#             on the Empirical Distribution Function.", Biometrika,
#             Vol. 66, Issue 3, Dec. 1979, pp 591-595.
_Avals_logistic = array([0.426, 0.563, 0.660, 0.769, 0.906, 1.010])
# From Richard A. Lockhart and Michael A. Stephens "Estimation and Tests of
#             Fit for the Three-Parameter Weibull Distribution"
#             Journal of the Royal Statistical Society.Series B(Methodological)
#             Vol. 56, No. 3 (1994), pp. 491-500, table 1. Keys are c*100
_Avals_weibull = [[0.292, 0.395, 0.467, 0.522, 0.617, 0.711, 0.836, 0.931],
                  [0.295, 0.399, 0.471, 0.527, 0.623, 0.719, 0.845, 0.941],
                  [0.298, 0.403, 0.476, 0.534, 0.631, 0.728, 0.856, 0.954],
                  [0.301, 0.408, 0.483, 0.541, 0.640, 0.738, 0.869, 0.969],
                  [0.305, 0.414, 0.490, 0.549, 0.650, 0.751, 0.885, 0.986],

Review comment (suggestion):
<img width="796" height="130" alt="image" src="https://github.com/user-attachments/assets/cc16ac6e-b6f3-404b-8103-5db90128c1ec" />


Code after:
    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)
    >>> 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.