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

Code before:
    #   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)
    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:
        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 * np.sqrt(df / ((1.0 - r + TINY)*(1.0 + r + TINY)))

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

        slope_stderr = np.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 * np.sqrt(ssxm + xmean**2)

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


Review comment (suggestion):
```suggestion
        prob = np.where(y[..., 0] == y[..., 1], 1.0, 0.0)
```

Code after:
        #        = 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,
        the result will broadcast correctly against the input array.
    nan_policy : {'propagate', 'omit', 'raise'}, default: 'propagate'
        Defines how to handle input NaNs.

        - ``propagate``: if a NaN is present in the axis slice (e.g. row) along
          which the statistic is computed, the corresponding entry of the output
          will be NaN.
        - ``omit``: NaNs will be omitted when performing the calculation.
          If insufficient data remains in the axis slice along which the
          statistic is computed, the corresponding entry of the output will be
          NaN.
        - ``raise``: if a NaN is present, a ``ValueError`` will be raised.

    dtype : dtype, optional
        Type to use in computing the mean. For integer inputs, the default is
        the default float type of the array library; for floating point inputs,
        the dtype is that of the input.

    Returns
    -------
    out : array
        The mean of each slice