322 lines
12 KiB
Python
322 lines
12 KiB
Python
import warnings
|
|
import numpy as np
|
|
from . import distributions
|
|
from .._lib._bunch import _make_tuple_bunch
|
|
from ._axis_nan_policy import _axis_nan_policy_factory
|
|
from ._stats_pythran import siegelslopes as siegelslopes_pythran
|
|
|
|
__all__ = ['_find_repeats', 'theilslopes', 'siegelslopes']
|
|
|
|
# This is not a namedtuple for backwards compatibility. See PR #12983
|
|
TheilslopesResult = _make_tuple_bunch('TheilslopesResult',
|
|
['slope', 'intercept',
|
|
'low_slope', 'high_slope'])
|
|
SiegelslopesResult = _make_tuple_bunch('SiegelslopesResult',
|
|
['slope', 'intercept'])
|
|
|
|
|
|
def _n_samples_optional_x(kwargs):
|
|
return 2 if kwargs.get('x', None) is not None else 1
|
|
|
|
|
|
@_axis_nan_policy_factory(TheilslopesResult, default_axis=None, n_outputs=4,
|
|
n_samples=_n_samples_optional_x,
|
|
result_to_tuple=lambda x, _: tuple(x), paired=True,
|
|
too_small=1)
|
|
def theilslopes(y, x=None, alpha=0.95, method='separate'):
|
|
r"""
|
|
Computes the Theil-Sen estimator for a set of points (x, y).
|
|
|
|
`theilslopes` implements a method for robust linear regression. It
|
|
computes the slope as the median of all slopes between paired values.
|
|
|
|
Parameters
|
|
----------
|
|
y : array_like
|
|
Dependent variable.
|
|
x : array_like or None, optional
|
|
Independent variable. If None, use ``arange(len(y))`` instead.
|
|
alpha : float, optional
|
|
Confidence degree between 0 and 1. Default is 95% confidence.
|
|
Note that `alpha` is symmetric around 0.5, i.e. both 0.1 and 0.9 are
|
|
interpreted as "find the 90% confidence interval".
|
|
method : {'joint', 'separate'}, optional
|
|
Method to be used for computing estimate for intercept.
|
|
Following methods are supported,
|
|
|
|
* 'joint': Uses np.median(y - slope * x) as intercept.
|
|
* 'separate': Uses np.median(y) - slope * np.median(x)
|
|
as intercept.
|
|
|
|
The default is 'separate'.
|
|
|
|
.. versionadded:: 1.8.0
|
|
|
|
Returns
|
|
-------
|
|
result : ``TheilslopesResult`` instance
|
|
The return value is an object with the following attributes:
|
|
|
|
slope : float
|
|
Theil slope.
|
|
intercept : float
|
|
Intercept of the Theil line.
|
|
low_slope : float
|
|
Lower bound of the confidence interval on `slope`.
|
|
high_slope : float
|
|
Upper bound of the confidence interval on `slope`.
|
|
|
|
See Also
|
|
--------
|
|
siegelslopes : a similar technique using repeated medians
|
|
|
|
Notes
|
|
-----
|
|
The implementation of `theilslopes` follows [1]_. The intercept is
|
|
not defined in [1]_, and here it is defined as ``median(y) -
|
|
slope*median(x)``, which is given in [3]_. Other definitions of
|
|
the intercept exist in the literature such as ``median(y - slope*x)``
|
|
in [4]_. The approach to compute the intercept can be determined by the
|
|
parameter ``method``. A confidence interval for the intercept is not
|
|
given as this question is not addressed in [1]_.
|
|
|
|
For compatibility with older versions of SciPy, the return value acts
|
|
like a ``namedtuple`` of length 4, with fields ``slope``, ``intercept``,
|
|
``low_slope``, and ``high_slope``, so one can continue to write::
|
|
|
|
slope, intercept, low_slope, high_slope = theilslopes(y, x)
|
|
|
|
References
|
|
----------
|
|
.. [1] P.K. Sen, "Estimates of the regression coefficient based on
|
|
Kendall's tau", J. Am. Stat. Assoc., Vol. 63, pp. 1379-1389, 1968.
|
|
.. [2] H. Theil, "A rank-invariant method of linear and polynomial
|
|
regression analysis I, II and III", Nederl. Akad. Wetensch., Proc.
|
|
53:, pp. 386-392, pp. 521-525, pp. 1397-1412, 1950.
|
|
.. [3] W.L. Conover, "Practical nonparametric statistics", 2nd ed.,
|
|
John Wiley and Sons, New York, pp. 493.
|
|
.. [4] https://en.wikipedia.org/wiki/Theil%E2%80%93Sen_estimator
|
|
|
|
Examples
|
|
--------
|
|
>>> import numpy as np
|
|
>>> from scipy import stats
|
|
>>> import matplotlib.pyplot as plt
|
|
|
|
>>> x = np.linspace(-5, 5, num=150)
|
|
>>> y = x + np.random.normal(size=x.size)
|
|
>>> y[11:15] += 10 # add outliers
|
|
>>> y[-5:] -= 7
|
|
|
|
Compute the slope, intercept and 90% confidence interval. For comparison,
|
|
also compute the least-squares fit with `linregress`:
|
|
|
|
>>> res = stats.theilslopes(y, x, 0.90, method='separate')
|
|
>>> lsq_res = stats.linregress(x, y)
|
|
|
|
Plot the results. The Theil-Sen regression line is shown in red, with the
|
|
dashed red lines illustrating the confidence interval of the slope (note
|
|
that the dashed red lines are not the confidence interval of the regression
|
|
as the confidence interval of the intercept is not included). The green
|
|
line shows the least-squares fit for comparison.
|
|
|
|
>>> fig = plt.figure()
|
|
>>> ax = fig.add_subplot(111)
|
|
>>> ax.plot(x, y, 'b.')
|
|
>>> ax.plot(x, res[1] + res[0] * x, 'r-')
|
|
>>> ax.plot(x, res[1] + res[2] * x, 'r--')
|
|
>>> ax.plot(x, res[1] + res[3] * x, 'r--')
|
|
>>> 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.")
|
|
# We copy both x and y so we can use _find_repeats.
|
|
y = np.array(y, dtype=float, copy=True).ravel()
|
|
if x is None:
|
|
x = np.arange(len(y), dtype=float)
|
|
else:
|
|
x = np.array(x, dtype=float, copy=True).ravel()
|
|
if len(x) != len(y):
|
|
raise ValueError("Array shapes are incompatible for broadcasting.")
|
|
if len(x) < 2:
|
|
raise ValueError("`x` and `y` must have length at least 2.")
|
|
|
|
# Compute sorted slopes only when deltax > 0
|
|
deltax = x[:, np.newaxis] - x
|
|
deltay = y[:, np.newaxis] - y
|
|
slopes = deltay[deltax > 0] / deltax[deltax > 0]
|
|
if not slopes.size:
|
|
msg = "All `x` coordinates are identical."
|
|
warnings.warn(msg, RuntimeWarning, stacklevel=2)
|
|
slopes.sort()
|
|
medslope = np.median(slopes)
|
|
if method == 'joint':
|
|
medinter = np.median(y - medslope * x)
|
|
else:
|
|
medinter = np.median(y) - medslope * np.median(x)
|
|
# Now compute confidence intervals
|
|
if alpha > 0.5:
|
|
alpha = 1. - alpha
|
|
|
|
z = distributions.norm.ppf(alpha / 2.)
|
|
# This implements (2.6) from Sen (1968)
|
|
_, nxreps = _find_repeats(x)
|
|
_, nyreps = _find_repeats(y)
|
|
nt = len(slopes) # N in Sen (1968)
|
|
ny = len(y) # n in Sen (1968)
|
|
# Equation 2.6 in Sen (1968):
|
|
sigsq = 1/18. * (ny * (ny-1) * (2*ny+5) -
|
|
sum(k * (k-1) * (2*k + 5) for k in nxreps) -
|
|
sum(k * (k-1) * (2*k + 5) for k in nyreps))
|
|
# Find the confidence interval indices in `slopes`
|
|
try:
|
|
sigma = np.sqrt(sigsq)
|
|
Ru = min(int(np.round((nt - z*sigma)/2.)), len(slopes)-1)
|
|
Rl = max(int(np.round((nt + z*sigma)/2.)) - 1, 0)
|
|
delta = slopes[[Rl, Ru]]
|
|
except (ValueError, IndexError):
|
|
delta = (np.nan, np.nan)
|
|
|
|
return TheilslopesResult(slope=medslope, intercept=medinter,
|
|
low_slope=delta[0], high_slope=delta[1])
|
|
|
|
|
|
def _find_repeats(arr):
|
|
# This function assumes it may clobber its input.
|
|
if len(arr) == 0:
|
|
return np.array(0, np.float64), np.array(0, np.intp)
|
|
|
|
# XXX This cast was previously needed for the Fortran implementation,
|
|
# should we ditch it?
|
|
arr = np.asarray(arr, np.float64).ravel()
|
|
arr.sort()
|
|
|
|
# Taken from NumPy 1.9's np.unique.
|
|
change = np.concatenate(([True], arr[1:] != arr[:-1]))
|
|
unique = arr[change]
|
|
change_idx = np.concatenate(np.nonzero(change) + ([arr.size],))
|
|
freq = np.diff(change_idx)
|
|
atleast2 = freq > 1
|
|
return unique[atleast2], freq[atleast2]
|
|
|
|
|
|
@_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="hierarchical"):
|
|
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.
|
|
|
|
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
|
|
--------
|
|
theilslopes : a similar technique without repeated medians
|
|
|
|
Notes
|
|
-----
|
|
With ``n = len(y)``, compute ``m_j`` as the median of
|
|
the slopes from the point ``(x[j], y[j])`` to all other `n-1` points.
|
|
``slope`` is then the median of all slopes ``m_j``.
|
|
Two ways are given to estimate the intercept in [1]_ which can be chosen
|
|
via the parameter ``method``.
|
|
The hierarchical approach uses the estimated slope ``slope``
|
|
and computes ``intercept`` as the median of ``y - slope*x``.
|
|
The other approach estimates the intercept separately as follows: for
|
|
each point ``(x[j], y[j])``, compute the intercepts of all the `n-1`
|
|
lines through the remaining points and take the median ``i_j``.
|
|
``intercept`` is the median of the ``i_j``.
|
|
|
|
The implementation computes `n` times the median of a vector of size `n`
|
|
which can be slow for large vectors. There are more efficient algorithms
|
|
(see [2]_) which are not implemented here.
|
|
|
|
For compatibility with older versions of SciPy, the return value acts
|
|
like a ``namedtuple`` of length 2, with fields ``slope`` and
|
|
``intercept``, so one can continue to write::
|
|
|
|
slope, intercept = siegelslopes(y, x)
|
|
|
|
References
|
|
----------
|
|
.. [1] A. Siegel, "Robust Regression Using Repeated Medians",
|
|
Biometrika, Vol. 69, pp. 242-244, 1982.
|
|
|
|
.. [2] A. Stein and M. Werman, "Finding the repeated median regression
|
|
line", Proceedings of the Third Annual ACM-SIAM Symposium on
|
|
Discrete Algorithms, pp. 409-413, 1992.
|
|
|
|
Examples
|
|
--------
|
|
>>> import numpy as np
|
|
>>> from scipy import stats
|
|
>>> import matplotlib.pyplot as plt
|
|
|
|
>>> x = np.linspace(-5, 5, num=150)
|
|
>>> y = x + np.random.normal(size=x.size)
|
|
>>> y[11:15] += 10 # add outliers
|
|
>>> y[-5:] -= 7
|
|
|
|
Compute the slope and intercept. For comparison, also compute the
|
|
least-squares fit with `linregress`:
|
|
|
|
>>> res = stats.siegelslopes(y, x)
|
|
>>> lsq_res = stats.linregress(x, y)
|
|
|
|
Plot the results. The Siegel regression line is shown in red. The green
|
|
line shows the least-squares fit for comparison.
|
|
|
|
>>> fig = plt.figure()
|
|
>>> ax = fig.add_subplot(111)
|
|
>>> ax.plot(x, y, 'b.')
|
|
>>> ax.plot(x, res[1] + res[0] * x, 'r-')
|
|
>>> ax.plot(x, lsq_res[1] + lsq_res[0] * x, 'g-')
|
|
>>> plt.show()
|
|
|
|
"""
|
|
if method not in ['hierarchical', 'separate']:
|
|
raise ValueError("method can only be 'hierarchical' or 'separate'")
|
|
y = np.asarray(y).ravel()
|
|
if x is None:
|
|
x = np.arange(len(y), dtype=float)
|
|
else:
|
|
x = np.asarray(x, dtype=float).ravel()
|
|
if len(x) != len(y):
|
|
raise ValueError("Array shapes are incompatible for broadcasting.")
|
|
if len(x) < 2:
|
|
raise ValueError("`x` and `y` must have length at least 2.")
|
|
|
|
dtype = np.result_type(x, y, np.float32) # use at least float32
|
|
y, x = y.astype(dtype), x.astype(dtype)
|
|
medslope, medinter = siegelslopes_pythran(y, x, method)
|
|
medslope, medinter = np.asarray(medslope)[()], np.asarray(medinter)[()]
|
|
return SiegelslopesResult(slope=medslope, intercept=medinter)
|