当前位置：首页 > article >正文

用矩阵和统计报告估计polynomial线性回归的系数python

article 2024/9/28 12:59:53

对于二次回归模型
y = beta0 + beta1 * x + beta2* x^2 + epsilon

在Python中，我们可以使用NumPy库来构建设计矩阵，并使用线性代数的方法来估计多项式线性回归的系数。

法⼀：矩阵

以下是一个简单的示例，展示了如何使用NumPy来估计一个二次多项式回归模型的系数：

import numpy as np
# 假设我们有一些数据点
X = np.array([1, 2, 3, 4, 5])  # 自变量
Y = np.array([2, 4, 5, 4, 5])  # 因变量
# 构建设计矩阵
# 对于二次多项式，我们需要添加X的平方项
X_poly = np.vstack([X**2, X, np.ones(len(X))]).T

# 使用NumPy的线性代数工具来计算系数
# np.linalg.inv() 计算逆矩阵，np.dot() 计算矩阵乘法
coefficients = np.linalg.inv(X_poly.T.dot(X_poly)).dot(X_poly.T).dot(Y)

# 打印系数
print("Coefficients:", coefficients)

在这个例子中，X_poly 是一个设计矩阵，它包含了 X 的平方项、X 的一次项和常数项。X_poly.T 是 X_poly 的转置。np.linalg.inv(X_poly.T.dot(X_poly)) 计算了 (X^T X)^(-1)，然后我们使用 .dot(X_poly.T).dot(Y) 来计算系数向量。

法二：LinearRegression

如果你想要使用更高级的库，比如 scikit-learn，你可以使用 PolynomialFeatures 来生成多项式特征，然后使用 LinearRegression 来拟合模型：

from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

def my_use_linear_regres(X, Y): 
    # 转换数据的形状
    X = X.reshape(-1, 1)

    # 创建一个PolynomialFeatures实例，degree指定多项式的度数
    poly = PolynomialFeatures(degree=2)   # change this!!!!!!!!!!!!!!!

    # 转换X
    X_poly = poly.fit_transform(X)

    # 创建线性回归模型实例
    model = LinearRegression()

    # 拟合模型
    model.fit(X_poly, Y)
    return model
model_mine = my_use_linear_regres(X, Y)
# 打印系数和截距
print("Coefficients:", model_mine.coef_)  # notice the first one is not meaningful!!! 
print("Intercept:", model_mine.intercept_)  # see intercept here!!!

在这个例子中，PolynomialFeatures 用于生成多项式特征，然后我们使用 LinearRegression 来拟合模型并提取系数。这种方法更加简洁，而且可以轻松地处理更复杂的多项式回归模型。

法三：report

在Python中，我们可以使用statsmodels库来执行线性回归分析，并获取一个详细的统计报告。

import numpy as np
import statsmodels.api as sm

# 假设我们有一些数据点
X = np.array([1, 2, 3, 4, 5])  # 自变量
Y = np.array([2, 4, 5, 4, 5])  # 因变量

def my_stat_report(X, Y):
        
    # 添加常数项以便拟合截距
    X = sm.add_constant(X)

    # 创建一个二次多项式模型的设计矩阵
    # 注意：这里我们直接使用X[:, 1]来表示x，X[:, 2]来表示x的平方
    X_poly_now = X[:, 1] ** 2
    # X_poly_now = np.column_stack((X[:, 1] ** 2, X[:, 1] ** 3))   # if there are more, use
    print(X_poly_now)
    # 将二次项和一次项添加到设计矩阵中
    X_design = np.column_stack((X, X_poly_now))
    print(X_design)
    # 拟合线性回归模型
    model = sm.OLS(Y, X_design).fit()
    return model

my_model = my_stat_report(X, Y)
print(my_model.summary())

说明
• Dep. Variable: 因变量的名称。
• R-squared: 决定系数，表示模型解释的变异性的比例。
• Adj. R-squared: 调整后的决定系数。
• F-statistic: F统计量，用于检验模型的整体显著性。
• Prob (F-statistic): F统计量的p值。
• Df Residuals: 残差的自由度。
• Df Model: 模型的自由度。
• Covariance Type: 协方差类型。
• coef: 回归系数。
• std err: 回归系数的标准误差。
• t: t统计量。
• P>|t|: t统计量的p值。
• [0.025 0.975]: 回归系数的95%置信区间。
这个统计报告提供了关于模型的详细统计信息，包括系数、标准误差、t统计量、p值和置信区间等。这些信息有助于评估模型的拟合质量和各个预测变量的显著性。

————

计算预测值 Y_pred 和残差平方和 (RSS)

# 计算预测值
Y_pred = model.predict(X_poly)

# 计算残差
residuals = Y - Y_pred

# 计算残差平方和 (RSS)
RSS = np.sum(residuals**2)

# 打印预测值和RSS
print("Predicted Y:", Y_pred)
print("Residual Sum of Squares (RSS):", RSS)

在这个代码中，model.predict(X_poly) 用于计算预测值。然后我们计算实际值 Y 和预测值 Y_pred 之间的残差，并将残差平方后求和得到 RSS。
RSS 越小，表示模型的拟合效果越好。RSS 是评估回归模型性能的一个常用指标。

RSE, MSE, TSS, R Square


# 计算总平方和 (TSS)
Y_mean = np.mean(Y)
TSS = np.sum((Y - Y_mean)**2)

# 计算均方误差 (MSE)
n = len(Y)  # 观测值数量
MSE = RSS / n 

# 计算剩余标准误差 (RSE)
N = len(Y)
p = 1
RSE = np.sqrt(RSS/(N-p-1))

# 计算 R-square
R_square = 1 - RSS / TSS

# 打印统计量
print("Residual Standard Error (RSE):", RSE)
print("Mean Squared Error (MSE):", MSE)
print("Total Sum of Squares (TSS):", TSS)
print("R-square:", R_square)