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

Code before:
    >>> ax.plot(x, lsq_res[1] + lsq_res[0] * x, 'g-')
    >>> plt.show()

    """
    if method not in ['joint', 'separate']:
        raise ValueError("method must be either 'joint' or 'separate'."
                         f"'{method}' is invalid.")

    y, x = xp_promote(y, x, force_floating=True, xp=np)
    x = np.arange(y.shape[-1], dtype=y.dtype) if x is None else x
    y, x = np.broadcast_arrays(y, x)

    if x.shape[-1] < 2 or y.shape[-1] < 2:  # only needed by test_axis_nan_policy
        raise ValueError("`x` and `y` must have length at least 2.")

    # Compute sorted slopes only when deltax > 0
    deltax = x[..., :, np.newaxis] - x[..., np.newaxis, :]
    deltay = y[..., :, np.newaxis] - y[..., np.newaxis, :]
    i = np.triu(np.ones(deltax.shape[-2:], dtype=bool), k=1)
    # with NumPy:
    deltax, deltay = deltax[..., i], deltay[..., i]
    # with array API, mask must be sole index, so we'll need to broadcast it.
    # indexing will ravel the array, so we'll need to reshape it after
    deltax[deltax == 0] = np.nan
    slopes = deltay / deltax
    slopes = np.sort(slopes, axis=-1)
    medslope = np.nanmedian(slopes, axis=-1)
    if method == 'joint':
        medinter = np.median(y - medslope * x, axis=-1)
    else:
        medinter = np.median(y, axis=-1) - medslope * np.median(x, axis=-1)
    # Now compute confidence intervals
    if alpha > 0.5:
        alpha = 1. - alpha

    z = stats.norm.ppf(alpha / 2.)
    # This implements (2.6) from Sen (1968)
    # we don't actually need ranks, so an enhancement could be to have
    # `rankdata` return only the second output. In the meantime, use the
    # least expensive `method`.
    _, nxreps = _rankdata(x, method='min', return_ties=True)
    _, nyreps = _rankdata(y, method='min', return_ties=True)
    nt = np.count_nonzero(np.isfinite(slopes), axis=-1, keepdims=True)  # N in Sen 1968
    ny = y.shape[-1]                                                    # n in Sen 1968
    # Equation 2.6 in Sen (1968):
    sigsq = 1/18. * (
        ny * (ny-1) * (2*ny+5)
        - np.sum(nxreps * (nxreps-1) * (2*nxreps + 5), axis=-1, keepdims=True)
        - np.sum(nyreps * (nyreps-1) * (2*nyreps + 5), axis=-1, keepdims=True))
    # Find the confidence interval indices in `slopes`
    sigma = np.sqrt(np.maximum(sigsq, 0.))

Review comment (suggestion):
This is only used for computing confidence intervals, which we don't do for `siegelslopes`, so it's moved until after returning the results for `siegelslopes`.

Code after:
    >>> ax.plot(x, res[1] + res[3] * x, 'r--')
    >>> ax.plot(x, lsq_res[1] + lsq_res[0] * x, 'g-')
    >>> plt.show()

    """
    return _robust_slopes(y, x=x, alpha=alpha, method=method, pfun='theilslopes')


@xp_capabilities(np_only=True)
@_axis_nan_policy_factory(SiegelslopesResult, default_axis=None, n_outputs=2,
                          n_samples=_n_samples_optional_x,
                          result_to_tuple=lambda x, _: tuple(x), paired=True,
                          too_small=1)
def siegelslopes(y, x=None, method='separate', *, axis=None):
    r"""
    Computes the Siegel estimator for a set of points (x, y).

    `siegelslopes` implements a method for robust linear regression
    using repeated medians (see [1]_) to fit a line to the points (x, y).
    The method is robust to outliers with an asymptotic breakdown point
    of 50%.

    Parameters
    ----------
    y : array_like
        Dependent variable.
    x : array_like or None, optional
        Independent variable. If None, use ``arange(len(y))`` instead.
    method : {'hierarchical', 'separate'}
        If 'hierarchical', estimate the intercept using the estimated
        slope ``slope`` (default option).
        If 'separate', estimate the intercept independent of the estimated
        slope. See Notes for details.
    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.

    Returns
    -------
    result : ``SiegelslopesResult`` instance
        The return value is an object with the following attributes:

        slope : float
            Estimate of the slope of the regression line.
        intercept : float
            Estimate of the intercept of the regression line.

    See Also
    --------