如何在Python中计算线性回归模型的AIC?

18

我想计算线性模型的AIC以比较它们的复杂性。 我是这样做的:

regr = linear_model.LinearRegression()
regr.fit(X, y)

aic_intercept_slope = aic(y, regr.coef_[0] * X.as_matrix() + regr.intercept_, k=1)

def aic(y, y_pred, k):
   resid = y - y_pred.ravel()
   sse = sum(resid ** 2)

   AIC = 2*k - 2*np.log(sse)

return AIC

但是我收到了一个除零错误

2个回答

21

sklearnLinearRegression用于预测很好,但像你发现的那样,它非常简单。(据说sklearn避开所有统计推断方面的东西。)

statsmodels.regression.linear_model.OLS具有属性AIC和其他一些预设属性。

然而,请注意,您需要手动向您的X矩阵添加一个单位向量,以在模型中包括截距。

from statsmodels.regression.linear_model import OLS
from statsmodels.tools import add_constant

regr = OLS(y, add_constant(X)).fit()
print(regr.aic)

如果您正在寻找一种另类的手写方式,同时仍然使用 sklearn,则可以在此处查看源代码


1

这里是一个AIC的示例实现,来自于之前回答中提供的链接

from statsmodels.regression.linear_model import OLS
from statsmodels.tools import add_constant
import numpy as np
from sklearn import linear_model

x = np.array([1,2,3,4,5])
y = np.array([0.3, 0.4, 0.4, 0.5, 0.6])

# 1 feature and constant
p = 1+1 

lr = linear_model.LinearRegression()
lr = lr.fit(x.reshape(-1, 1), y)
pr = lr.predict(x.reshape(-1, 1))

def llf_(y, X, pr):
    # return maximized log likelihood
    nobs = float(X.shape[0])
    nobs2 = nobs / 2.0
    nobs = float(nobs)
    resid = y - pr
    ssr = np.sum((resid)**2)
    llf = -nobs2*np.log(2*np.pi) - nobs2*np.log(ssr / nobs) - nobs2
    return llf

def aic(y, X, pr, p):
    # return aic metric
    llf = llf_(y, X, pr)
    return -2*llf+2*p

regr = OLS(y, add_constant(x.reshape(-1, 1))).fit()
print(regr.aic)
print(aic(y, x, pr, p))

输出:

-18.903519181693923 # OLS AIC
-18.903519181693916 # our AIC

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接