scipy/scipy — pull request #3574 by tomflannaghan (2014-04-24 16:16:00+00:00)
Title: Fixes multiple bugs in mstats.theilslopes

This fixes multiple bugs in `scipy.stats.mstats.theilslopes`. I believe that the original code is implementing the method given in Sen (1968) for the slope and confidence intervals (this is the standard method), and I have therefore corrected the implementation to follow Sen (1968) faithfully. I have also added a test that highlights some of the failures of the original code.

An outline of the changes:

1. The original implementation erroneously always masked the last element of the input data (line 699 in the original code), and so the last element was always omitted from the remaining calculations.

2. The code retained the masked values in the `slopes` array, which caused problems when the confidence intervals are calculated by array index. I now use `x.compressed()` and `y.compressed()` to remove masked values.

3. Slopes were being computed for repeated values of x, not in accordance with the Sen (1968) method. This is fixed by replacing the routine for computing `slopes` with one that only includes values where `deltax > 0`.

4. The original calculation of `sigsq` does not divide the second term by 18 (line 722 in original code). It also includes a `yties` term for handling repeated `y` values (line 723 in original code) that does not appear in Sen (1968) equation (2.6).

5. The confidence intervals are computed by computing indices in the `slopes` array. In Sen (1968), these indices are 1-indexed, but python is 0-indexed. Therefore, we must subtract 1 from the values computed by Sen (1968) to index the `slopes` array correctly (see lines 728-729 in the new code).