Странно как-то вы вычисляете доверительный интервал для коэффициентов. Не знаю, на каких учебниках вы основываетесь, но вот посмотрите тут:

http://mcimeer.narod.ru/data/t5/t5_2.html

Она отличается от той, которую вы пытаетесь применить и несколько отличается от той, которую (для коэффициента b) привел уважаемый Pak Uula

Если к штатным возможностям вы относите scipy и numpy, то доверительные интервалы для линейной регрессии можно найти несколькими способами.

Какой бы способ вы ни выбрали, для вычисления доверительного интервала нужно знать сам параметр (например, a) и его среднеквадратичное отклонение (пусть будет a_err).

Доверительный интервал с уровнем доверия alpha вычисляют по распределению Стьюдента:

conf_int = scipy.stats.t.interval(1-alpha, df=n-2, loc=a, scale=a_err)

Если нужно найти только полуширину интервала - ту величину, которая ставится после знака ±, то она вычисляется так:

plus_minus = abs(sps.t.ppf(alpha/2, n-2))*a_err

Теперь как найти параметры и их ошибки.

Средствами linregress

scipy.stats.linregress - специализированный метод вычисления линейной регрессии.

import scipy.stats as sps

n = len(x)
lin_model = sps.linregress(x, y)
a,b = lin_model.slope, lin_model.intercept
# оценка ср.кв. ошибки для a и b
a_err, b_err = lin_model.stderr, lin_model.intercept_stderr
# Доверительный интервал для alpha=5%
a_conf = sps.t.interval(0.95, df = n-2, loc=a, scale=a_err)
b_conf = sps.t.interval(0.95, df = n-2, loc=b, scale=b_err)

print(f"a = {a:0.4f}, α=5% [{a_conf[0]:0.4f} - {a_conf[1]:0.4f}]")
print(f"b = {b:0.4f}, α=5% [{b_conf[0]:0.4f} - {b_conf[1]:0.4f}]")

Результат для 100 точек на прямой y=0.5x+2 со случайной ошибкой sigma=0.5:

a = 0.5291, α=5% [0.4971 - 0.5610]
b = 1.8568, α=5% [1.6718 - 2.0418]

Универсальный инструмент curve_fit

В scipy есть универсальный инструмент приближения набора точек заданной моделью scipy.optimize.curve_fit. Эта функция ищет наилучший набор параметров методом наименьших квадратов для любого вида моделей, не только линейных. Помимо оптимальных значений параметров функция возвращает ковариационную матрицу, диагональные элементы которой дают оценку дисперсии параметров.

import numpy as np
import scipy.stats as sps
import scipy.optimize as spo

def linear(x, a,b):
    return a*x+b

((a,b), cov) = spo.curve_fit(linear, xdata=x, ydata=y)
a_err, b_err = np.sqrt(np.diag(cov))
a_conf = sps.t.interval(0.95, df = n-2, loc=a, scale=a_err)
b_conf = sps.t.interval(0.95, df = n-2, loc=b, scale=b_err)

Результат для того же набора данных:

a = 0.5291, α=5% [0.4971 - 0.5610]
b = 1.8568, α=5% [1.6718 - 2.0418]

Как видно, результат тот же, что для специализированного средства.

Прямое вычисление

Можно напрямую вычислить параметры линейной регрессии по формулам

# Вычисление параметров модели
sum_x = x.sum()
sum_y = y.sum()
sum_xy = (x*y).sum()
sum_x_sq = (x*x).sum()
a = (n*sum_xy - sum_x*sum_y)/(n*sum_x_sq - sum_x*sum_x)
b = (sum_y*sum_x_sq - sum_x*sum_xy)/(n*sum_x_sq - sum_x*sum_x)

# вычисление ошибки параметров
u = y - (a*x+b)
u_avg = np.mean(u)
sigma_square = 1.0/(n-2)*np.sum((u - u_avg)**2)
x_mean = np.mean(x)
dx_square = np.sum((x-x_mean)**2)

a_err = np.sqrt(sigma_square/dx_square)
b_err = np.sqrt(sigma_square*(1.0/n + np.mean(x)**2/dx_square))

Замеры производительности

Линейная регрессия для 100 точек

1. Прямое вычисление: 51.6 µs ± 2.47
2. linregress : 207 µs ± 2.67 µs
3. curve_fit: 237 µs ± 8.29 µs

Результат линейной регрессии на графике.

Пример

Пример загружен в github как jupyter notebook.