pute the seasonal slopes as:

>>> from scipy import stats
>>> intra_slope, inter_slope = stats.mstats.sen_seasonal_slopes(x)

If we define a function to compute all slopes between observations within
a season:

>>> def dijk(yi):
...     n = len(yi)
...     x = np.arange(n)
...     dy = yi - yi[:, np.newaxis]
...     dx = x - x[:, np.newaxis]
...     # we only want unique pairs of distinct indices
...     mask = np.triu(np.ones((n, n), dtype=bool), k=1)
...     return dy[mask]/dx[mask]

then element ``i`` of ``intra_slope`` is the median of ``dijk[x[:, i]]``:

>>> i = 2
>>> np.allclose(np.median(dijk(x[:, i])), intra_slope[i])
True

and ``inter_slope`` is the median of the values returned by ``dijk`` for
all seasons:

>>> all_slopes = np.concatenate([dijk(x[:, i]) for i in range(x.shape[1])])
>>> np.allclose(np.median(all_slopes), inter_slope)
True

Because the data are randomly generated, we would expect the median slopes
to be nearly zero both within and across all seasons, and indeed they are:

>>> intra_slope.data
array([ 0.00124504, -0.00277761, -0.00221245, -0.00036338])
>>> inter_slope
-0.0010511779872922058

TFra