对一组数据进行拟合构建模型的过程实际上是从一个model family中找到最优参数的过程。书中从一个简单的线性模型开始带领我们进入model的世界，提供的练习题帮助我们理解记忆model的知识，此处仅做一个简单的解答。
问题1：线性模型的一个缺点是对异常值敏感，因为root-mean-squared deviation中包含了一个平方项。用线性模型拟合下面的模拟数据，并多次生成模拟数据观察model的变化。
# Q1:生成模拟数据集sim1a(如下)，拟合一个线性模型并可视化结果。重复运行多次，注意model的结果有何不同。
# 生成模拟数据集
# rt(n,df,ncp):随机生成一组数量为n,自由度为df的t分布的数据集
# ncp: non-centrality parameter delta, 如果忽略，则使用中心t分布
# 构建模型,同时也计算了model的预测值
model1 <- function(coefs,data){
coefs[1] + coefs[2] * data$x
}
# 计算预测值与实际值之间的距离来评价模型的优劣，评价方法：均方根误差(root-mean-squared deviation)
measure_distance <- function(coefs,data){
diff <- data$y - model1(coefs,data)
sqrt(mean(diff^2))
}
# ggplot2可视化(使用root-mean-squared)
visual1_ggplot2 <- function(model,data,text1="lm()"){
deviation <- measure_distance(model,data)
fit_formula <- paste0("y=",round(model[1],2),"+",round(model[2],2),"*x(root-mean-squared=",round(deviation),")")
p<- ggplot(data,aes(x=x,y=y)) + geom_point() + geom_abline(intercept=model[1],slope=model[2])
p<- p + labs(title=fit_formula,subtitle = text1) + theme(plot.title = element_text(hjust = 0.5))
return(p)
}
# 可视化方式一：
#将绘图的默认设置存储起来
op<-par(no.readonly = T)
par(mfrow=c(3,3))
for(i in 1:9){
sim1a <- tibble(
x = rep(1:10, each = 3),
y = x * 1.5 + 6 + rt(length(x), df = 2)
)
sim1a_mod <- lm(y~x,data=sim1a)
deviation <- measure_distance(c(sim1a_mod$coefficients[1],sim1a_mod$coefficients[2]),sim1a)
fit_formula <- paste0("y=",round(sim1a_mod$coefficients[1]),"+",round(sim1a_mod$coefficients[2]),"*x(mse=",round(deviation),")")
plot(sim1a$x,sim1a$y,pch=19,xlab="x",ylab="y",main=fit_formula)
abline(a=sim1a_mod$coefficients[1],b=sim1a_mod$coefficients[2])
}
# 恢复默认设置
par(op)
# 可视化方式二：
library(ggplot2)
library(ggpubr)
pictures <- vector(mode="list",length = 9)
for(i in 1:9){
sim1a <- tibble(
x = rep(1:10, each = 3),
y = x * 1.5 + 6 + rt(length(x), df = 2)
)
mod <- lm(y~x,data=sim1a)
pictures[[i]] <- visual1_ggplot2(mod$coefficients,sim1a)
}
ggarrange(plotlist=pictures,ncol=3,nrow=3)
Q1结果
有异常值的数据，拟合效果较差？
问题2：使线性模型更稳健的方式使使用一个更稳健的距离度量，例如使用平均绝对距离(mean-absolute distance)，使用optim()拟合模型并比较与线性模型的差别。
######################################################################################################################################
# Q2:使线性模型更稳健的方式使使用一个更稳健的距离度量，例如使用平均绝对距离(mean-absolute distance)
# 使用optim()拟合模型并比较与线性模型的差别。
# 使用平均绝对距离评价模型优劣
measure_distance1 <- function(coefs,data){
diff <- data$y - model1(coefs,data)
mean(abs(diff))
}
# ggplot2可视化(使用mean-absolute)
visual2_ggplot2 <- function(model,data,text1="lm()"){
deviation <- measure_distance1(model,data)
fit_formula <- paste0("y=",round(model[1],2),"+",round(model[2],2),"*x(mean-absolute=",round(deviation),")")
p<- ggplot(data,aes(x=x,y=y)) + geom_point() + geom_abline(intercept=model[1],slope=model[2])
p<- p + labs(title=fit_formula,subtitle = text1) + theme(plot.title = element_text(hjust = 0.5))
return(p)
}
# 生成模拟数据
sim1 <- tibble(
x = rep(1:10, each = 3),
y = x * 1.5 + 6 + rt(length(x), df = 2)
)
# 使用optim()方法，距离度量为root-mean-squared deviation
best1 <- optim(c(0, 0), measure_distance, data = sim1)
p1<- visual1_ggplot2(best1$par,sim1,"optim()")
# 使用optim()方法，距离度量为mean-absolute
best2 <- optim(c(0, 0), measure_distance1, data = sim1)
p2<- visual2_ggplot2(best2$par,sim1,"optim()")
# 使用lm()方法，距离度量为root-mean-squared deviation
lm1 <- lm(y~x,data=sim1)
p3<- visual1_ggplot2(lm1$coefficients,sim1,"lm()")
# 使用lm()方法，距离度量为mean-absolute
lm2 <- lm(y~x,data=sim1)
p4<- visual2_ggplot2(lm2$coefficients,sim1,"lm()")
ggarrange(p1,p2,p3,p4,ncol=2,nrow = 2)
Q2的结果
使用mean-absolute距离度量，optim()方法比lm()方法的拟合效果更好？
问题3::使用数值优化的一个挑战是只能保证找到一个局部最优。优化如下的三参数模型有什么问题吗？
######################################################################################################################################
# Q3:使用数值优化的一个挑战是只能保证找到一个局部最优。优化如下的三参数模型有什么问题吗？
model1 <- function(a, data) {
a[1] + data$x * a[2] + a[3]
}
# 本问题奇奇怪怪，可可爱爱。
不能理解问题3中a[1]+a[3]不就是一个常数吗？为什么会有3个参数？
网友评论