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

Code before:
    >>> print(f"R-squared: {res.rvalue**2:.6f}")
    R-squared: 0.717533

    Plot the data along with the fitted line:

    >>> plt.plot(x, y, 'o', label='original data')
    >>> plt.plot(x, res.intercept + res.slope*x, 'r', label='fitted line')
    >>> plt.legend()
    >>> plt.show()

    Calculate 95% confidence interval on slope and intercept:

    >>> # Two-sided inverse Students t-distribution
    >>> # 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

    """
    xp = array_namespace(x, y)  # needed for vecdot pre NumPy 2.0
    x, y = xp_promote(x, y, force_floating=True, xp=xp)

    TINY = 1.0e-20

    # _axis_nan_policy decorator ensures that `axis=-1`
    n = x.shape[-1]
    xmean = xp.mean(x, axis=-1, keepdims=True)
    ymean = xp.mean(y, axis=-1, 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=-1, xp=xp)
    y_ = _demean(y, ymean, axis=-1, xp=xp, precision_warning=False)
    xmean = xp.squeeze(xmean, axis=-1)
    ymean = xp.squeeze(ymean, axis=-1)

    ssxm = xp.vecdot(x_, x_, axis=-1) / n
    ssym = xp.vecdot(y_, y_, axis=-1) / n
    ssxym = xp.vecdot(x_, y_, axis=-1) / n

    # R-value
    #   r = ssxym / sqrt( ssxm * ssym )
    degenerate = (ssxm == 0.0) | (ssym == 0.0)
    NaN = xp.asarray(xp.nan, dtype=ssxym.dtype)

Review comment (suggestion):
```suggestion
    xp = array_namespace(x, y)
```

Code after:

    >>> print(f"R-squared: {res.rvalue**2:.6f}")
    R-squared: 0.717533

    Plot the data along with the fitted line:

    >>> plt.plot(x, y, 'o', label='original data')
    >>> plt.plot(x, res.intercept + res.slope*x, 'r', label='fitted line')
    >>> plt.legend()
    >>> plt.show()

    Calculate 95% confidence interval on slope and intercept:

    >>> # Two-sided inverse Students t-distribution
    >>> # 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

    """
    xp = array_namespace(x, y)
    x, y = xp_promote(x, y, force_floating=True, xp=xp)

    TINY = 1.0e-20

    # _axis_nan_policy decorator ensures that `axis=-1`
    n = x.shape[-1]
    xmean = xp.mean(x, axis=-1, keepdims=True)
    ymean = xp.mean(y, axis=-1, 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=-1, xp=xp)
    y_ = _demean(y, ymean, axis=-1, xp=xp, precision_warning=False)
    xmean = xp.squeeze(xmean, axis=-1)
    ymean = xp.squeeze(ymean, axis=-1)

    ssxm = xp.vecdot(x_, x_, axis=-1) / n
    ssym = xp.vecdot(y_, y_, axis=-1) / n
    ssxym = xp.vecdot(x_, y_, axis=-1) / n

    # R-value
    #   r = ssxym / sqrt( ssxm * ssym )
    degenerate = (ssxm == 0.0) | (ssym == 0.0)