The seasonal generalization of Sen's slope computes the slopes between all pairs of values within a "season" (column) of a 2D array. It returns an array containing the median of these "within-season" slopes for each season (the Theil-Sen slope estimator of each season), and it returns the median of the within-season slopes across all seasons (the seasonal Kendall slope estimator).

The slopes d_{ijk} within season i are:

d_{ijk} = \frac{x_{ij} - x_{ik}}
                    {j - k}

for pairs of distinct integer indices j, k of x.

Element i of the returned intra_slope array is the median of the d_{ijk} over all j < k; this is the Theil-Sen slope estimator of season i. The returned inter_slope value, better known as the seasonal Kendall slope estimator, is the median of the d_{ijk} over all i, j, k.

Computes seasonal Theil-Sen and Kendall slope estimators.

import numpy as np
rng = np.random.default_rng()
x = rng.random(size=(100, 4))

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

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]

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

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

intra_slope.data

inter_slope