《统计方法与机器学习》实践作业2¶
10211900416 郭夏辉
D-Day：Oct. 18¶
背景描述¶
汽车发动机在测功机上产生的制动马力被认为是发动机转速(每分钟转数，rpm)、燃料的道路辛烷值和发动机压缩值的函数，我们在实验室里进行实验，研究它们的函数关系。
数据描述¶
| 变量名 | 变量含义 | 变量类型 | 变量取值范围 |
|---|---|---|---|
| rpm | 发动机转速 | 连续变量 | $R^+$ |
| Road_Octane_Number | 道路辛烷值 | 连续变量 | $R+$ |
| Compression | 压缩值 | 连续变量 | $R^+$ |
| Brake_Horsepower | 制动马力 | 连续变量 | $R^+$ |
任务¶
注：这里使用 $\alpha=0.05$ 的显著性水平。
- 请用多元线性回归模型，描述制动马力和发动机转速、道路辛烷值以及压缩值之间的函数关系。
- 分别将数据中心化、标准化之后，比较参数估计的异同，并进行评述（提示：可以结合理论课的课件）。
- 从模型显著性、参数显著性以及残差分析三个角度，分析多元线性回归模型是否合理。
- 若取发动机转速为3000转/min，道路辛烷值为90，发动机压缩值为100时，分别给出制动马力值的置信区间和预测区间。
import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt
from statsmodels.formula.api import ols
from scipy.stats import f
from scipy.stats import t
from sklearn import datasets, linear_model
from sklearn.metrics import mean_squared_error, r2_score
from sklearn import preprocessing
alpha = 0.05
p = 3
n = 12
df = pd.read_csv("Project_3.csv")
col = {
'Brake Horsepower':'y',
'rpm':'xi1',
'Road Octane Number':'xi2',
'Compression':'xi3'
}
df.rename(columns=col,inplace=True)
df['intercept'] = 1 # 截距要乘的数，就是X矩阵第一列，为1
df = df[['intercept','xi1', 'xi2', 'xi3','y']]
data = df.values
X = data[:, 0:(p+1)]
y = data[:, -1]
print(df)
#print(X)
#print(y)
intercept xi1 xi2 xi3 y 0 1 2000 90 100 225 1 1 1800 94 95 212 2 1 2400 88 110 229 3 1 1900 91 96 222 4 1 1600 86 100 219 5 1 2500 96 110 278 6 1 3000 94 98 246 7 1 3200 90 100 237 8 1 2800 88 105 233 9 1 3400 86 97 224 10 1 1800 90 100 223 11 1 2500 89 104 230
1.请用多元线性回归模型，描述制动马力和发动机转速、道路辛烷值以及压缩值之间的函数关系。¶
线性回归模型：$$y=\beta_{0}+\beta_{1}x_{1}+...+\beta_{p}x_{p}+\epsilon$$
线性回归模型的数据化写法：$$y_{i}=\beta_{0}+\beta_{1}x_{i1}+...+\beta_{p}x_{ip}+\epsilon_{i}$$
可以这样来写 $$\mathbf{y} =\mathbf{X} \mathbf{\beta} + \mathbf{\epsilon} $$
自变量$$\mathbf{X} = \begin{pmatrix} 1 & x_{11} & x_{12} & x_{13}\\ 1 & x_{21} & x_{22} & x_{23}\\ \vdots & \vdots & \vdots & \vdots\\ 1 & x_{n1} & x_{n2} & x_{n3}\\ \end{pmatrix}$$
响应变量$$\mathbf{y} = \begin{pmatrix}y_1\\y_2\\\vdots\\ y_n\end{pmatrix}$$
待估参数$$\mathbf{\beta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \end{pmatrix}$$
随机误差$$\mathbf{\epsilon} = \begin{pmatrix}\epsilon_1\\\epsilon_2\\\vdots\\ \epsilon_n\end{pmatrix}$$
方法一、直接代入公式¶
根据课程所学知识，$\beta$的最小二乘估计是：
$$\hat{\beta} = (X^T X)^{-1}X^T y$$
然后，拟合值定义为：
$$\hat{y} = X(X^T X)^{-1}X^T y$$
beta_hat_1 = np.linalg.inv(X.T @ X) @ (X.T @ y)
print('线性回归方程直接代公式为:')
print(f"y_hat = {beta_hat_1[0]} + {beta_hat_1[1]} * x1 + {beta_hat_1[2]} * x2 + {beta_hat_1[3]} * x3")
线性回归方程直接代公式为: y_hat = -266.03121172235114 + 0.010713207851841267 * x1 + 3.1348062584375143 * x2 + 1.8674094338437044 * x3
方法二、调用statsmodels包中的方法¶
# statsmodels
model1=ols('y~xi1 + xi2 + xi3',df).fit()
beta_hat_2 = model1.params
print(f"参数估计值: {beta_hat_2}")
参数估计值: Intercept -266.031212 xi1 0.010713 xi2 3.134806 xi3 1.867409 dtype: float64
方法三、调用scikit-learn包中的方法¶
model2 = linear_model.LinearRegression()
X_without_intercept = X[:,1:(p+1)]
model2.fit(X_without_intercept, y)
beta_hat_3 = np.append(np.array(model2.intercept_),model2.coef_)
print(f"参数估计值: {beta_hat_3}")
参数估计值: [-2.66031212e+02 1.07132079e-02 3.13480626e+00 1.86740943e+00]
可以看到，这三种方法得到的结果是一样的。
2.分别将数据中心化、标准化之后，比较参数估计的异同，并进行评述（提示：可以结合理论课的课件）。¶
2.1 中心化¶
中心化的含义：
将矩阵$\mathbf{X}$中的每一个$x_{ij}$和向量$\mathbf{y}$中的每一个$y_{i}$都对该列的均值作差，一定注意是按列计算，而不是按行计算！我最开始学习时候就有点混淆。
$$x_{ij}^*=x_{ij}-\frac{1}{n}\sum_{i=1}^{n} x_{ij}$$ $$y_{i}^*=y_{i}-\frac{1}{n}\sum_{i=1}^{n} y_{i}$$
X_mean = []
for k in range(p + 1):
X_mean.append(np.mean(data[:, k])) # 第k列的均值
y_mean = np.mean(data[:, -1])
X_cent = X - X_mean
y_cent = y - y_mean
df1 = pd.DataFrame(X_cent,columns=['intercept', 'xi1', 'xi2','xi3'])
df1['y'] = y_cent
model3 = ols('y ~ xi1 + xi2 + xi3', df1).fit()
beta_centhat_1 = model3.params
print(f"中心化后，参数估计值: {beta_centhat_1}")
中心化后，参数估计值: Intercept -4.440892e-16 xi1 1.071321e-02 xi2 3.134806e+00 xi3 1.867409e+00 dtype: float64
与题目1的结果对比，我发现除了$\beta_{0}$的估计变为0(虽然输出的并不是0，但是10的-16次方这个数量级可以认为是趋于0了)，其他的参数的估计并没有变。如果要说为什么会出现这样的情况，还要回归理论课程，这里就不赘述了。
但我想补充一下为何会出现这样的情况之理解，在中心化之后，线性模型中斜率项没有变，只是改变了截距项，只消除了截距对响应变量带来的影响不会导致除了$\beta_{0}$之外的变化。
然后还有一种对数据进行中心化的方法，就不用我在那先按列求方差再做差了——直接调用sklearn的preprocessing方法。
X_cent1 = preprocessing.scale(X, with_mean = True, with_std=False)
y_cent1 = preprocessing.scale(y, with_mean = True, with_std=False)
print(X_cent)
print(X_cent1)
print(y_cent)
print(y_cent1)
[[ 0.00000000e+00 -4.08333333e+02 -1.66666667e-01 -1.25000000e+00] [ 0.00000000e+00 -6.08333333e+02 3.83333333e+00 -6.25000000e+00] [ 0.00000000e+00 -8.33333333e+00 -2.16666667e+00 8.75000000e+00] [ 0.00000000e+00 -5.08333333e+02 8.33333333e-01 -5.25000000e+00] [ 0.00000000e+00 -8.08333333e+02 -4.16666667e+00 -1.25000000e+00] [ 0.00000000e+00 9.16666667e+01 5.83333333e+00 8.75000000e+00] [ 0.00000000e+00 5.91666667e+02 3.83333333e+00 -3.25000000e+00] [ 0.00000000e+00 7.91666667e+02 -1.66666667e-01 -1.25000000e+00] [ 0.00000000e+00 3.91666667e+02 -2.16666667e+00 3.75000000e+00] [ 0.00000000e+00 9.91666667e+02 -4.16666667e+00 -4.25000000e+00] [ 0.00000000e+00 -6.08333333e+02 -1.66666667e-01 -1.25000000e+00] [ 0.00000000e+00 9.16666667e+01 -1.16666667e+00 2.75000000e+00]] [[ 0.00000000e+00 -4.08333333e+02 -1.66666667e-01 -1.25000000e+00] [ 0.00000000e+00 -6.08333333e+02 3.83333333e+00 -6.25000000e+00] [ 0.00000000e+00 -8.33333333e+00 -2.16666667e+00 8.75000000e+00] [ 0.00000000e+00 -5.08333333e+02 8.33333333e-01 -5.25000000e+00] [ 0.00000000e+00 -8.08333333e+02 -4.16666667e+00 -1.25000000e+00] [ 0.00000000e+00 9.16666667e+01 5.83333333e+00 8.75000000e+00] [ 0.00000000e+00 5.91666667e+02 3.83333333e+00 -3.25000000e+00] [ 0.00000000e+00 7.91666667e+02 -1.66666667e-01 -1.25000000e+00] [ 0.00000000e+00 3.91666667e+02 -2.16666667e+00 3.75000000e+00] [ 0.00000000e+00 9.91666667e+02 -4.16666667e+00 -4.25000000e+00] [ 0.00000000e+00 -6.08333333e+02 -1.66666667e-01 -1.25000000e+00] [ 0.00000000e+00 9.16666667e+01 -1.16666667e+00 2.75000000e+00]] [ -6.5 -19.5 -2.5 -9.5 -12.5 46.5 14.5 5.5 1.5 -7.5 -8.5 -1.5] [ -6.5 -19.5 -2.5 -9.5 -12.5 46.5 14.5 5.5 1.5 -7.5 -8.5 -1.5]
可以发现是一样的，sklearn的一些函数真的简洁又好用！
2.2 标准化¶
标准化处理数据：
$$x_{ij}^{**}=\frac{x_{ij}^*}{\sqrt{\sum_{i=1}^n (x_{ij}-\bar{x_{j}})^2}}$$
$$y_{i}^{**}=\frac{y_{i}^*}{\sqrt{\sum_{i=1}^n (y_{i}-\bar{y})^2}}$$
注意这里还是像中心化一样一列一列地去计算。
这次我为了简便，直接调用preprocessing方法来标准化吧。
X_scale = preprocessing.scale(X, with_mean = True, with_std=True)
y_scale = preprocessing.scale(y, with_mean = True, with_std=True)
df2 = pd.DataFrame(X_scale,columns=['intercept', 'xi1', 'xi2','xi3'])
df2['y'] = y_scale
model4 = ols('y ~ xi1 + xi2 + xi3', df2).fit()
beta_scalehat_1 = model4.params
print(f"标准化后，参数估计值: {beta_scalehat_1}")
标准化后，参数估计值: Intercept -2.775558e-17 xi1 3.757094e-01 xi2 5.793325e-01 xi3 5.477351e-01 dtype: float64
与题目1的结果对比，我发现：
$\beta_{0}$的估计变为0(虽然输出的并不是0，但是10的-17次方这个数量级可以认为是趋于0了)
$\beta_{1}$的估计值变大了，$\beta_{2}和\beta_{3}$的估计值变小了。
但是在原来$\beta_{2} \ge \beta_{3} \ge \beta_{1}$
现在中心化之后还是$\beta_{2} \ge \beta_{3} \ge \beta_{1}$
标准化消除了变量的量纲对响应变量产生的影响，但并不会影响到自变量本身大小对响应变量的影响。
3.从模型显著性、参数显著性以及残差分析三个角度，分析多元线性回归模型是否合理。¶
3.1 模型显著性¶
要检验模型的显著性，即检验： $$H_0: \beta_{1}=\beta_{2}=\beta_{3}=0 v.s. H_1: 存在一个\beta_{i}\neq0 i=1,2,3$$
如果$H_0$成立，说明马力和发动机转速、燃料的道路辛烷值和发动机压缩值并没有什么关系，模型不合适。
检验模型的显著性使用$F$检验。检验统计量：
$$F_{0}=\frac{SS_R/p}{SS_E/(n-p-1)}$$
其中 $$SS_R=\sum_{i=1}^n (\hat{y_{i}}-\bar{y})^2$$ $$SS_E=\sum_{i=1}^n (y_{i}-\hat{y_{i}})^2$$
根据理论课程的内容，当$H_0$成立时 $F_{0} \sim F(p,n-p-1)$
在显著性水平$\alpha$下，拒绝域为 $\{F_0 > F_{1-\alpha}(p,n-p-1)\}$，$p$值为$p=P(F\geqslant F_0)$
y_hat_1 = model1.fittedvalues
SSE = sum((y - y_hat_1 ) ** 2)
SST = sum((y - y_mean ) ** 2)
SSR = sum((y_hat_1 - y_mean) ** 2)
sigma2 = SSE / (n-p-1)
F0 = (SSR / p) / (SSE / (n - p - 1))
print(F0)
11.1159636363607
难道我真的还要费劲地去算$SS_R$和$SS_E$再求得$F_0$嘛？不是的。我直接可以通过模型的报告信息就能看到相关的取值。
model1.summary()
C:\learnAI\anaconda3\lib\site-packages\scipy\stats\_stats_py.py:1736: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=12
warnings.warn("kurtosistest only valid for n>=20 ... continuing "
| Dep. Variable: | y | R-squared: | 0.807 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.734 |
| Method: | Least Squares | F-statistic: | 11.12 |
| Date: | Fri, 13 Oct 2023 | Prob (F-statistic): | 0.00317 |
| Time: | 14:45:30 | Log-Likelihood: | -40.708 |
| No. Observations: | 12 | AIC: | 89.42 |
| Df Residuals: | 8 | BIC: | 91.36 |
| Df Model: | 3 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | -266.0312 | 92.674 | -2.871 | 0.021 | -479.737 | -52.325 |
| xi1 | 0.0107 | 0.004 | 2.390 | 0.044 | 0.000 | 0.021 |
| xi2 | 3.1348 | 0.844 | 3.712 | 0.006 | 1.188 | 5.082 |
| xi3 | 1.8674 | 0.535 | 3.494 | 0.008 | 0.635 | 3.100 |
| Omnibus: | 0.392 | Durbin-Watson: | 1.043 |
|---|---|---|---|
| Prob(Omnibus): | 0.822 | Jarque-Bera (JB): | 0.230 |
| Skew: | -0.282 | Prob(JB): | 0.891 |
| Kurtosis: | 2.625 | Cond. No. | 9.03e+04 |
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 9.03e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
在上表中可以很清楚地看到F检验统计量为11.12，p值0.00317，小于显著性水平，故我们拒绝原假设，认为这个模型是显著的。
3.2 参数显著性¶
在多元线性回归中，回归方程显著并不意味着每个自变量对因变量的影响都显著。我们希望从回归方程中剔除那些次要的、可有可无的自变量，建立更为简化的回归方程。因此我们需要对每个自变量进行显著性检验。如果某个$x_j$对$y$的影响不显著，那么在回归模型中，其对应的回归系数设置为0.
$$H_{0j}: \beta_j=0 v.s. H_{1j}:\beta_j\neq 0, j\in\{1,2..p\} $$ 如果原假设成立，我们就认为自变量$x_j$对响应变量没有影响，即影响不显著。
我们一般用$t$检验法来检验上述假设检验问题。
检验统计量 $t_j=\frac{\hat{\beta_j}}{\sqrt{c_{jj} \hat{\sigma}}}$
${\hat{\sigma}}^2 = \frac{1}{n-p-1}SS_E$
由于$H_0$成立时$t_j \sim t(n-p-1)$显著性水平$\alpha$下拒绝域为$\{|t_j|\geqslant t_{1-\alpha/2}(n-p-1)\}$
C = np.linalg.inv(np.dot(X.T, X))
t0 = []
for i in range(p + 1):
t0.append(beta_hat_1[i] / (np.sqrt(C[i][i] * sigma2))) # 求t值
print(t0)
print('t的临界值为：',t.ppf(1 - alpha / 2, n - p - 1))
[-2.8706239247834833, 2.3896037316292653, 3.7123120444773643, 3.4935816668577373] t的临界值为： 2.3060041350333704
其实在模型显著性检验时，通过那张大表格我们就能看到三个变量的t检验统计量及p值了。它们的p值分别为0.044，0.006，0.008，都小于显著性水平，都落入了拒绝域中，各回归系数都满足显著性。
3.3 残差分析¶
这个其实就是第一个实验内容的一部分。我们要验证这个模型的正态假设是否成立，即验证 $$\epsilon_i\sim N(0,\sigma^2),i\in [1,n]\bigcap Z$$且各$\epsilon_i$相互独立。
我们要检验残差的独立性、方差齐性和正态性。
# 计算残差
df['res'] = df["y"]-(beta_hat_1[0]+beta_hat_1[1]*df["xi1"]+beta_hat_1[2]*df["xi2"]+beta_hat_1[3]*df["xi3"])
print(df['res'] )
0 0.731289 1 -13.328247 2 -11.958476 3 3.137442 4 11.555798 5 10.891754 6 2.213675 7 -0.124560 8 -2.906712 9 2.874252 10 0.873931 11 -3.960146 Name: res, dtype: float64
3.3.1 独立性¶
# Durbin-Watson检验
def durbin_watson(residuals):
nume = sum(np.diff(residuals.T) ** 2)
deno = sum(residuals ** 2)
return nume / deno
dw = durbin_watson(df['res'])
print('Durbin-Watson检验的检验统计量为：',dw)
Durbin-Watson检验的检验统计量为： 1.0431120863346202
根据实验一，一种较为粗略的判断方法：$DW$一般介于0到4之间，如果$DW$值接近于2，我们一般可以认为残差序列不具有相关性。
$DW$求出来是1.04，在0到4之间。虽然距离2有一点距离，但是可以近似地认为残差之间相互独立
3.3.2 方差齐性¶
本来我想着用Levene检验来检验方差齐性呢，但是看到只有一组的残差数据，就知道这个没必要做了。
3.3.3 正态性¶
# Shapiro-Wilk检验
SW, pVal1 = stats.shapiro(df['res'])
print(SW)
print(pVal1)
if pVal1 > alpha:
print('接受H0')
else:
print('p-value < alpha, 落入拒绝域，拒绝H0')
0.9290957450866699 0.3706260919570923 接受H0
根据实际结果，统计量SW为 0.9291，接近 1 且P值为 0.371，大于指定的显著性水平 0.05。 故认为残差来自服从正态分布的总体。
3.4 决定系数的检验¶
虽然这个实验并没有要求来做这个，但是在课程中讲过一些关于决定系数的知识。
定义样本决定系数为 $R^{2}=\frac{S S_{R}}{S S_{T}}=1-\frac{S S_{E}}{S S_{T}}$
$R^{2}$ 的取值在 [0,1] 区间内。
$R^{2}$ 越接近 $1,$ 表明回归拟合的效果越好。
$R^{2}$ 越接近 $0,$ 表明回归拟合的效果越差。
然后我想尝试着做一下。
R2 = SSR / SST
print('决定系数R2=', R2)
决定系数R2= 0.8065197565314398
虽然这个决定系数距离1还有些距离，但是0.8整体来说还不错，说明模型的回归拟合效果还行。
4.若取发动机转速为3000转/min，道路辛烷值为90，发动机压缩值为100时，分别给出制动马力值的置信区间和预测区间。¶
在线性回归模型中，响应变量的点预测值就等于它在该点的估计值
$$\hat{y}=\mathbf{x^T}(\mathbf{X^TX})^{-1}\mathbf{X^Ty}$$
注意这里的$\mathbf{X}$矩阵和$\mathbf{x}$向量别搞混了。
$E[y]$置信区间为$\hat{y} \pm t_{1-\alpha/2}(n-p-1) \hat{\sigma} \sqrt{\mathbf{x}^T (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{x}}$
$y$预测区间为$\hat{y} \pm t_{1-\alpha/2}(n-p-1) \hat{\sigma} \sqrt{1+\mathbf{x}^T (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{x}}$
其中$\hat{\sigma}^2 = (n-p-1) \sum_{i=1}^n e_i^2 = \frac{1}{n-p-1}SS_E$
sigma2 = SSE / (n-p-1)
tval = t.ppf(1 - alpha / 2, n - p - 1)
x0 = [1,3000,90,100]
x0 = np.array(x0)
y0hat = np.dot(x0.T, beta_hat_1)
delta = tval * np.sqrt(sigma2) * np.sqrt(x0.T @ C @ x0)
delta1 = tval * np.sqrt(sigma2) * np.sqrt(1+x0.T @ C @ x0)
Y = [y0hat - delta, y0hat + delta]
Y1 = [y0hat - delta1, y0hat + delta1]
print("点估计=",y0hat)
print("E[y]置信区间=",Y)
print("y预测区间=",Y1)
点估计= 234.9819184769194 E[y]置信区间= [226.24574501813427, 243.7180919357045] y预测区间= [212.8622464639742, 257.1015904898646]
不要忽略两者计算公式不同！计算得到的$y$的预测区间要宽于$E[y]$的置信区间。