主页 / 问题 / 1283898
Accepted
Юрий
Asked:2022-05-20 02:18:53 +0000 UTC2022-05-20 02:18:53 +0000 UTC

Python中回归模型参数的置信区间

  • 772

有2个样本X，Y，需要建立一个模型Y=aX+b+eps。有必要找到参数 a,b 的置信区间。

根据学生的相应公式，我自己编写了区间搜索。（teta-参数矩阵，252-样本量，假设正态分布，即2个参数，可靠性水平a=0.05）

a_left = teta[0]-stats.t.ppf(0.975,250)*math.sqrt(eps2.sum()*C[0,0]/250)
a_right = teta[0]+stats.t.ppf(0.975,250)*math.sqrt(eps2.sum()*C[0,0]/250)
b_left = teta[1]-stats.t.ppf(0.975,250)*math.sqrt(eps2.sum()*C[1,1]/250)
b_right = teta[1]+stats.t.ppf(0.975,250)*math.sqrt(eps2.sum()*C[1,1]/250)

但是，我需要（违背我的意愿）只使用库的标准功能（stats、sklearn、numpy、scipy 等）。这些库有这样的内置功能吗？

目前只能找到模型中的参数

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn import metrics

model = LinearRegression().fit(x, y)
y_pred = model.predict(x)
print("a,b= ",  model.coef_[0], model.intercept_)
python

2 个回答

  1. Best Answer

    如果您将和归因于常规可能性，则可以通过多种方式找到线性回归的置信区间。scipynumpy

    无论您选择哪种方法，要计算置信区间，您都需要知道参数本身（例如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

    图上线性回归的结果。

    线性回归

    例子

    该示例作为 jupyter notebook 上传到 github。

    • 3

  2. 奇怪的是，您以某种方式计算系数的置信区间。我不知道你是基于什么教程，但看看这里：

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

    它与您尝试应用的不同，并且与受人尊敬的 Pak Uula 给出的（系数 b）有所不同

    • 0

相关问题