机器学习—线性回归与R模拟

1、数据模拟

# 数据模拟
set.seed(1)
N <- 10000
x1 <- 10+(20-10)*runif(N) #生成10000个10-20的随机数
x2 <- 10+(50-10)*runif(N) #生成10000个10-50的随机数
e <- rnorm(N,0,2) #生成1000个正态分布均值0，方差2
y <- 0.3*x1+ 0.4*x2 - 1.2*x1*x2 - 0.6*I(x1^2)+ e
summary(lm(y~x1*x2+ I(x1^2)))

# Call:
#   lm(formula = y ~ x1 * x2 + I(x1^2))
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -8.6198 -1.3567 -0.0348  1.3293  7.4732 
# 
# Coefficients:
#             Estimate Std. Error   t value Pr(>|t|)    
# (Intercept) -0.1747179  0.6391113    -0.273 0.784570    
# x1           0.2978905  0.0812568     3.666 0.000248 ***
# x2           0.4040485  0.0089673    45.058  < 2e-16 ***
# I(x1^2)     -0.5994130  0.0026345  -227.525  < 2e-16 ***
# x1:x2       -1.2001345  0.0005887 -2038.484  < 2e-16 ***
#   ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.982 on 9995 degrees of freedom
# Multiple R-squared:  0.9999,	Adjusted R-squared:  0.9999 
# F-statistic: 4.305e+07 on 4 and 9995 DF,  p-value: < 2.2e-16

mydata <- data.frame(y = y,x1 = x1,x2 = x2) #合并为数据框

其中模型
$$y=0.3x_1+0.4x_2-1.2x_1{x}_2-0.6x_1^2+e$$
为数据生成过程（DGP）。通过OLS回归得到
$$y=-0.17-00.2978x_1+0.4040x_2-1.200x_1{x}_2-0.5994x_1^2+e$$
为总体回归模型。

---

2、验证集法

验证集法将样本数据一分为二，将其中一部分作为为训练集（train sample），一部分作为测试集(test sample)；用训练集估计模型的参数，并以估计后的模型预测测试集，并计算验证集误差（validation set error），以此估计测试集误差(test error)

# 验证集法：从10000个数据随机分为两部分70%为训练集，30%为测试集
set.seed(2)
train <- sample(N,0.7*N)  # 随机抽取70%
fit <- lm(y~x1*x2+ I(x1^2),subset = train,data = mydata)
summary(fit)
# Call:
# lm(formula = y ~ x1 * x2 + I(x1^2), data = mydata, subset = train)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -6.7103 -1.3641 -0.0341  1.3063  7.4762 
# 
# Coefficients:
#              Estimate Std. Error   t value Pr(>|t|)    
# (Intercept) -0.1736355  0.7631043    -0.228  0.82001    
# x1           0.3089380  0.0969150     3.188  0.00144 ** 
# x2           0.3996129  0.0106803    37.416  < 2e-16 ***
# I(x1^2)     -0.6000714  0.0031394  -191.141  < 2e-16 ***
# x1:x2       -1.1998799  0.0007038 -1704.854  < 2e-16 ***
#   ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.985 on 6995 degrees of freedom
# Multiple R-squared:  0.9999,	Adjusted R-squared:  0.9999 
# F-statistic: 2.988e+07 on 4 and 6995 DF,  p-value: < 2.2e-16

predict_train_data <- predict(fit,mydata[-train,]) # 测试集预测值
y.test <- mydata$y[-train] # 测试集实际值
MSE <- mean((predict_train_data- y.test)^2) # 计算均方误差
print(MSE)
# 3.898763

将7000样本观测回归，预测3000样本的y值，计算均方误差为3.8987，说明模型外严外推能力较强，但验证集法局限是：将70%观测作为训练集，损失了30%的样本观测，模型估计效果存在信息损失。一种方法是通过设定不同随机种子，再不同随机种子下计算平均MSE。但整体而言，验证集法与实际结果存在偏差。在ML中，大量使用K折交叉验证法。

---

3 、K折交叉验证法

K折交叉验证法是将全样本划分为K个子集，一般$K=5,10$。再每次保留一个子集作为测试集，其余$K-1$个子集作为训练集，共计算$K$个均方误差(MSE）。设第$k$折的均方误差为$MS{E}_k$，则$K$折交叉验证误差为（平均MSE）
$$C{V}_k=\frac{1}{K}\sum _{k=1}^{K}MS{E}_k$$
更极端的，当折数$K=n$($n$为样本容量)时，即将每一个样本观测作为测试集，其余$n-1$个观测作为训练集，计算$n$次交叉验证误差，这种方法称为留一交叉验证法。

library(boot)
fit <- glm(y~x1*x2+ I(x1^2),data = mydata)
set.seed(1)
cv.err <- cv.glm(mydata,fit,K = 10) #K为样本划分折数，通常取K = 5或10
names(cv.err)
# "call"  "K"  "delta" "seed" 
cv.err$delta  # delta为均方误差；结果有两个，第一个是普通MSE，第二十校正后的MSE
# 3.930578 3.930255   

# 设定不同随机种子
set.seed(10)
cv.err <- cv.glm(mydata,fit,K = 10)
cv.err$delta
# 3.929462 3.929198

# 留一交叉验证法
# 将每一个样本观测作为测试集，另外9999作为训练集，
# 并预测测试集合，计算均方误差MSE;需要等待一些时间
cv.err <- cv.glm(mydata,fit)
cv.err$delta
# 3.928419 3.928419

从上述几种模型评估方法看，验证集法、K折交叉验证法和留一交叉验证法的MSE差别不大且较小，可以判断线性回归模型具有良好的外推或预测能力。

参考文献

陈强.机器学习及R应用[M].高等教育出版社

---

-END-