Why do the two sides of theilslopes() not meet neatly at the median value? I originally posted this question to StackExchange and was requested to move it here, as it relates to domain knowledge rather than code writing.
I’ve been playing with the theilslopes() function in Python to determine if the trends in two datasets are significantly different. While exploring this question I thought it would be useful to plot the probability distributions of the two datasets. I did this by running 50 different alpha levels through the function and plotting the upper and lower bounds. After doing this I was struck by the gap in the middle of the plot where the trend suddenly jumps from one side to the other through the median value.
Why are the upper and lower confidence levels not equal to the median value for alpha = 0.50?

n_pts = 50
rng = np.random.default_rng()
x1 = rng.random(n_pts)
y1 = (np.random.rand()-0.5)*x1 + rng.random(n_pts)
alphas = [x / 100 for x in range(1, 51)]
temp = alphas.copy()
temp.reverse()
alphas_all = np.array(alphas + temp[1:])
slopes = np.full([len(alphas) * 2 - 1], np.nan)

for ii, alpha in enumerate(alphas):
    a1 = fun_c_theilslopes(y1, x1, alpha=alpha)
    slopes[ii] = a1[2]
    slopes[-1 - ii] = a1[3]
a1 = fun_c_theilslopes(y1, x1, alpha=0.05)

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(4.5, 2.0), sharex=False, sharey=False)
axes[0].plot(x1, y1, "ok", ms=1, label="original data")
axes[0].plot(x1, a1[1] + (a1[0] * x1), "r")
intercept = a1[1] + (np.median(x1) * (a1[0] - a1[3]))
axes[0].plot([min(x1), max(x1)], intercept + (a1[3] * np.array([min(x1), max(x1)])),
             c='r', linestyle='--', dashes=(4, 2), label='upper')
intercept = a1[1] + (np.median(x1) * (a1[0] - a1[2]))
axes[0].plot([min(x1), max(x1)], intercept + (a1[2] * np.array([min(x1), max(x1)])),
             c='r', linestyle='--', dashes=(4, 2), label='lower')
axes[0].set_ylabel("y-value")
axes[0].set_xlabel("x-value")
axes[0].set_title("Original data points and slope with \n alpha = 0.05 confidence limits")

axes[1].plot(slopes, alphas_all, "ob", ms=1, label="original data")
axes[1].plot(slopes[[4, -5]], alphas_all[[4, -5]], "or", ms=3, label="original data")
axes[1].plot([a1[0], a1[0]], [0, 0.5],
             c="r", linestyle="-", label="median")
axes[1].set_ylabel("alpha")
axes[1].set_xlabel("slope")
axes[1].set_title("Upper and lower confidence \n limits for each alpha level")
plt.pause(5)
plt.show(block=True)
plt.close()

Following @Glen_b's answer (the gap is the 50% CI) I have plotted slopes corresponding to confidence intervals down to 1%. The theilslopes() function treats alpha = 0.05 and alpha = 0.95 the same - both are used to calculate the 95% confidence interval.


A symmetric 95% interval would have (1-0.95)/2 = 2.5% in each tail (outside the interval) and 95% within the interval. A symmetric 50% interval would have 25% in each tail and 50% within the interval (discreteness issues aside, which might change things slightly).
So, naturally, you should tend to expect them to have a gap.
