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

Code before:
    >>> # p - probability, df - degrees of freedom
    >>> from scipy.stats import t
    >>> tinv = lambda p, df: abs(t.ppf(p/2, df))

    >>> ts = tinv(0.05, len(x)-2)
    >>> print(f"slope (95%): {res.slope:.6f} +/- {ts*res.stderr:.6f}")
    slope (95%): 1.453392 +/- 0.743465
    >>> print(f"intercept (95%): {res.intercept:.6f}"
    ...       f" +/- {ts*res.intercept_stderr:.6f}")
    intercept (95%): 0.616950 +/- 0.544475

    """
    np = array_namespace(x, y)  # needed for vecdot pre NumPy 2.0

    TINY = 1.0e-20

    axis = -1
    n = x.shape[axis]
    xmean = np.mean(x, axis=axis, keepdims=True)
    ymean = np.mean(y, axis=axis, keepdims=True)

    # Average sums of square differences from the mean
    #   ssxm = mean( (x-mean(x))^2 )
    #   ssxym = mean( (x-mean(x)) * (y-mean(y)) )
    x_ = _demean(x, xmean, axis=axis, xp=np)
    y_ = _demean(y, ymean, axis=axis, xp=np, precision_warning=False)
    xmean = np.squeeze(xmean, axis=-1)
    ymean = np.squeeze(ymean, axis=-1)

    ssxm = np.vecdot(x_, x_, axis=axis) / n
    ssym = np.vecdot(y_, y_, axis=axis) / n
    ssxym = np.vecdot(x_, y_, axis=axis) / n

    # R-value
    #   r = ssxym / sqrt( ssxm * ssym )
    degenerate = (ssxm == 0.0) | (ssym == 0.0)
    r = xpx.apply_where(
        ~degenerate,
        (ssxym, ssxm, ssym),
        lambda ssxym, ssxm, ssym: np.clip(ssxym / np.sqrt(ssxm * ssym), -1.0, 1.0),
        lambda ssxym, ssxm, ssym: np.where(ssxym==0, np.nan, 0.0)
    )

    slope = ssxym / ssxm
    intercept = ymean - slope*xmean
    if n == 2:
        # handle case when only two points are passed in
        prob = np.where(y[..., 0] == y[..., 1], 1.0, 0.0)
        slope_stderr = np.zeros_like(r)
        intercept_stderr = np.zeros_like(r)
    else:

Review comment (suggestion):
```suggestion
    x_ = _demean(x, xmean, axis=-1, xp=np)
    y_ = _demean(y, ymean, axis=-1, xp=np, precision_warning=False)
```

Code after:
    if n == 2:
        # handle case when only two points are passed in
        one = xp.asarray(1.0, dtype=r.dtype)
        prob = xp.where(y[..., 0] == y[..., 1], one, 0.0)
        slope_stderr = xp.zeros_like(r)
        intercept_stderr = xp.zeros_like(r)
    else:
        df = n - 2  # Number of degrees of freedom
        # n-2 degrees of freedom because 2 has been used up
        # to estimate the mean and standard deviation
        t = r * xp.sqrt(df / ((1.0 - r + TINY)*(1.0 + r + TINY)))

        dist = _SimpleStudentT(xp.asarray(df, dtype=t.dtype))
        prob = _get_pvalue(t, dist, alternative, xp=xp)
        prob = prob[()] if prob.ndim == 0 else prob

        slope_stderr = xp.sqrt((1 - r**2) * ssym / ssxm / df)

        # Also calculate the standard error of the intercept
        # The following relationship is used:
        #   ssxm = mean( (x-mean(x))^2 )
        #        = ssx - sx*sx
        #        = mean( x^2 ) - mean(x)^2
        intercept_stderr = slope_stderr * xp.sqrt(ssxm + xmean**2)

    return LinregressResult(slope=slope[()], intercept=intercept[()], rvalue=r[()],
                            pvalue=prob[()], stderr=slope_stderr[()],
                            intercept_stderr=intercept_stderr[()])


def _xp_mean(x, /, *, axis=None, weights=None, keepdims=False, nan_policy='propagate',
             dtype=None, xp=None):
    r"""Compute the arithmetic mean along the specified axis.

    Parameters
    ----------
    x : real array
        Array containing real numbers whose mean is desired.
    axis : int or tuple of ints, default: None
        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.
    weights : real array, optional
        If specified, an array of weights associated with the values in `x`;
        otherwise ``1``. If `weights` and `x` do not have the same shape, the
        arrays will be broadcasted before performing the calculation. See
        Notes for details.
    keepdims : boolean, optional
        If this is set to ``True``, the axes which are reduced are left
        in the result as dimensions with length one. With this option,