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

Code before:
    # order necessary conditions are satisfied. If not, raise an error.
    # Here, use `m` for the shape parameter to be consistent with [7]
    # and avoid confusion with `c` as defined in [7].
    n = len(x)
    m, u, s = params

    def dnllf_dm(m, u):
        # Partial w.r.t. shape w/ optimal scale. See [7] Equation 5.
        xu = x-u
        return (1/m - (xu**m*np.log(xu)).sum()/(xu**m).sum()
                + np.log(xu).sum()/n)

    def dnllf_du(m, u):
        # Partial w.r.t. loc w/ optimal scale. See [7] Equation 6.
        xu = x-u
        return (m-1)/m*(xu**-1).sum() - n*(xu**(m-1)).sum()/(xu**m).sum()

    def get_scale(m, u):
        # Partial w.r.t. scale solved in terms of shape and location.
        # See [7] Equation 7.
        return ((x-u)**m/n).sum()**(1/m)

    def dnllf(params):
        # Partial derivatives of the NLLF w.r.t. parameters, i.e.
        # first order necessary conditions for MLE fit.
        return [dnllf_dm(*params), dnllf_du(*params)]

    suggestion = ("Maximum likelihood estimation is known to be challenging "
                  "for the three-parameter Weibull distribution. Consider "
                  "performing a custom goodness-of-fit test using "
                  "`scipy.stats.monte_carlo_test`.")

    if np.allclose(u, np.min(x)) or m < 1:
        # The critical values provided by [7] don't seem to control the
        # Type I error rate in this case. Error out.
        message = ("Maximum likelihood estimation has converged to "
                   "a solution in which the location is equal to the minimum "
                   "of the data, the shape parameter is less than 2, or both. "
                   "The table of critical values in [7] does not "
                   "include this case. " + suggestion)
        raise ValueError(message)

    try:
        # Refine the MLE / verify that first-order necessary conditions are
        # satisfied. If so, the critical values provided in [7] seem reliable.
        with np.errstate(over='raise', invalid='raise'):
            res = optimize.root(dnllf, params[:-1])

        message = ("Solution of MLE first-order conditions failed: "
                   f"{res.message}. `anderson` cannot continue. " + suggestion)
        if not res.success:

Review comment (suggestion):
<img width="887" height="81" alt="image" src="https://github.com/user-attachments/assets/45423eb2-df70-4378-b1f7-eaa5bde9a67e" />


Code after:
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)
    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 = xpx.at(a)[..., -2].set(anm1)
        a = xpx.at(a)[..., 1].set(-anm1)
    a = xpx.at(a)[..., -1].set(an)
    a = xpx.at(a)[..., -0].set(-an)
    W = _swilk_w(y, a, xp=xp)

    # 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 = math.exp(1.3822 - 0.77857*u + 0.062767*u**2 - 0.0020322*u**3)
        gW = -xp.log(gamma - xp.log(1 - W))
    else:
        u = math.log(n)
        mu = -1.5861 - 0.31082*u - 0.083751*u**2 + 0.0038915*u**3
        sigma = math.exp(-0.4803 - 0.082676*u + 0.0030302*u**2)
        gW = xp.log(1 - W)

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