R - 核支持向量机预测模型构建详解：以短波辐射通量shortw
2019-12-30 本文已影响0人
TroyShen
目录
- 0.问题导入
- 1.示例数据
- 2.导入示例数据
- 3.将示例数据按照70%/30%拆分为训练组与验证组
- 4.构建核支持向量机预测模型（SVR）
- 5.模型训练结果与实测数据相关性及RMSE分析（图1）
- 6.模型训练结果与实测数据时间路径分析（图2）
- 7.SVR模型预测结果与验证集实测相关性及RMSE分析（图3）
- 8.SVR模型预测结果与验证集实测数据时间路径分析（图4）
- 9.模型模拟残差分析(图5)
- 10.总结
- 11.本篇所使用的软件包（没有的需要通过install.packages('包名')进行安装）
- 12.致谢
0. 问题导入
今天我们通过建立短波辐射通量与温度及蒸散发的关系来详细阐述核支持向量机预测模型（SVR）在面向地理及时间序列数据建模方面的应用。
1.示例数据
本篇数据与R - 广义加性模型（GAM）构建详解：以归一化植被指数NDVI为例这篇文章数据一致，时间长度为2002-01 ~ 2015-09，时间分辨率为月。
点我下载示例数据
示例数据预览：
head(data)
ndvi soil1 soil2 soil3 soil4 rain
1 0.237984169 -0.01210715 0.03579731 0.1269299 0.07318894 -0.01543584
2 0.370455335 0.38147139 0.31089661 0.2241396 0.10204067 0.20701857
3 0.331657733 0.41044975 0.48385978 0.4471074 0.25112199 0.62105802
4 0.216662956 0.32583872 0.41198999 0.4231082 0.42613716 0.37216417
5 0.054132382 0.24177292 0.20540345 0.2979310 0.43549429 0.06553887
6 -0.005636952 0.41268755 0.29207486 0.2508858 0.37087816 0.25502620
longwave shortwave root_sm evap temper gpp
1 0.059414987 0.215758745 0.06890271 -0.07747205 0.009909431 -0.04072053
2 0.009142641 0.244385277 0.31129426 0.23793998 0.172678808 0.18329118
3 -0.097150000 0.353491078 0.48706357 0.59985033 0.314583437 0.24478460
4 -0.031527285 0.355970841 0.42948289 0.51469995 0.348457057 0.47172457
5 -0.291598633 0.297255464 0.27643746 0.38420614 0.326270291 0.37802032
6 0.010729154 -0.009709685 0.32028271 0.31684306 0.119881673 0.15986512
2. 导入示例数据
setwd('/Users/jerseyshen/Documents/JianShu_Project/20191230')
data = read.csv("data.csv",
header = T)
data = data[,-1]
date = seq(as.Date('2002-01-01'),
as.Date('2015-12-01'),
'1 month')[1:165]
3. 将示例数据按照70%/30%拆分为训练组与验证组
train = data[1:115,]
test = data[116:165,]
date_train = date[1:115]
date_test = date[116:165]
4. 构建核支持向量机预测模型（SVR）
本例采用的SVR模型中，回归模型为：eps-regression ； kennel为：RBF(Radial Base Function)。
model_init = svm(shortwave ~ temper+evap,
data = train,
kernel = 'radial')
5. 模型训练结果与实测数据相关性及RMSE分析（图1）
train_simu = as.numeric(model_init$fitted)
pl_point_df = data.frame(SIMU = as.numeric(train_simu),
TRAIN = train$shortwave)
label_df = data.frame(
x = c(-0.5,-0.5,-0.5),
y = c(0.3,0.25,0.2),
label = c(paste0("COR = ", round(cor(train_simu,train$shortwave),2)),
paste0("p = ", formatC(cor.test(train_simu,train$shortwave)$p.value)),
paste0("RMSE = ", round(rmse(train_simu,train$shortwave),3)))
)
p1_cor = ggplot()+
geom_point(data = pl_point_df,
aes(x = SIMU, y= TRAIN),
size = 5,
color = 'blue',
shape = 1)+
geom_abline(intercept = 0,slope = 1,size = 1)+
geom_text(data = label_df,aes(x = x,y = y,label = label),size = 6,color = 'black',
hjust = 0)+
theme_bw()+
theme(
axis.text = element_text(face = 'bold',color = 'black',size = 12),
axis.title = element_text(face = 'bold',color = 'black',size = 14, hjust = .5)
)+
xlab('SIMULATION RESULT')+
ylab('REAL SHORTWAVE')
png('plot1.png',
height = 20,
width = 20,
units = 'cm',
res = 800)
print(p1_cor)
dev.off()
图1 模型训练结果与实测数据相关性及RMSE分析
6. 模型训练结果与实测数据时间路径分析（图2）
pl_line_df = data.frame(
Date = date_train,
SIMU = train_simu,
REAL = train$shortwave
)
pl_line_df = melt(pl_line_df,'Date')
p1_line = ggplot()+
geom_line(data = pl_line_df,
aes(x = Date, y= value,color= variable,linetype = variable),
size= 1)+
scale_color_manual(values = c('green','blue'))+
scale_linetype_manual(values = c('solid','dashed'))+
theme_bw()+
theme(
axis.text = element_text(face = 'bold',color = 'black',size = 12),
axis.title = element_text(face = 'bold',color = 'black',size = 14, hjust = .5),
legend.text = element_text(face = 'bold',color = 'black',size = 12),
legend.title = element_blank(),
legend.position = 'bottom',
legend.direction = 'horizontal'
)+
xlab('Date')+
ylab('Indices')
png('plot2.png',
height = 10,
width = 20,
units = 'cm',
res = 800)
print(p1_line)
dev.off()
图2 模型训练结果与实测数据时间路径分析
7. SVR模型预测结果与验证集实测相关性及RMSE分析（图3）
test_simu = predict(model_init,
test)
pl_point_df2 = data.frame(SIMU = as.numeric(test_simu),
TEST = test$shortwave)
label_df2 = data.frame(
x = c(-0.5,-0.5,-0.5),
y = c(0.3,0.25,0.2),
label = c(paste0("COR = ", round(cor(test_simu,test$shortwave),2)),
paste0("p = ", formatC(cor.test(test_simu,test$shortwave)$p.value)),
paste0("RMSE = ", round(rmse(test_simu,test$shortwave),3)))
)
########
p2_cor = ggplot()+
geom_point(data = pl_point_df2,
aes(x = SIMU, y= TEST),
size = 5,
color = 'blue',
shape = 1)+
geom_abline(intercept = 0,slope = 1,size = 1)+
geom_text(data = label_df2,aes(x = x,y = y,label = label),size = 6,color = 'black',
hjust = 0)+
theme_bw()+
theme(
axis.text = element_text(face = 'bold',color = 'black',size = 12),
axis.title = element_text(face = 'bold',color = 'black',size = 14, hjust = .5)
)+
xlab('SIMULATION RESULT')+
ylab('REAL SHORTWAVE')
png('plot3.png',
height = 20,
width = 20,
units = 'cm',
res = 800)
print(p2_cor)
dev.off()
图3 SVR模型预测结果与验证集实测相关性及RMSE分析
8. SVR模型预测结果与验证集实测数据时间路径分析（图4）
pl_line_df2 = data.frame(
Date = date_test,
SIMU = test_simu,
REAL = test$shortwave
)
pl_line_df2 = melt(pl_line_df2,'Date')
p1_line2 = ggplot()+
geom_line(data = pl_line_df2,
aes(x = Date, y= value,color= variable,linetype = variable),
size= 1)+
scale_color_manual(values = c('green','blue'))+
scale_linetype_manual(values = c('solid','dashed'))+
theme_bw()+
theme(
axis.text = element_text(face = 'bold',color = 'black',size = 12),
axis.title = element_text(face = 'bold',color = 'black',size = 14, hjust = .5),
legend.text = element_text(face = 'bold',color = 'black',size = 12),
legend.title = element_blank(),
legend.position = 'bottom',
legend.direction = 'horizontal'
)+
xlab('Date')+
ylab('Indices')
png('plot4.png',
height = 10,
width = 20,
units = 'cm',
res = 800)
print(p1_line2)
dev.off()
图4 SVR模型预测结果与验证集实测数据时间路径分析
9. 模型模拟残差分析(图5)
由图5可见，SVR模型在训练集与验证集的模拟残差大致服从于均值为0, 标准差为0.1 的正态分布，证明模型模拟与预测能力强。
TRAIN_RES = model_init$fitted - train$shortwave
TEST_RES = test_simu - test$shortwave
df_res1 = data.frame(RES = TRAIN_RES, Type = "TRAIN_RES")
df_res2 = data.frame(RES = TEST_RES, Type = "TEST_RES")
df_res = rbind(df_res1,df_res2)
label_df_res = data.frame(
x = c(0.3,0.3,0.3,0.3),
y = c(500,450,400,350),
labels = c(
paste0('TRAIN_MEAN = ',round(mean(TRAIN_RES),2)),
paste0('TRAIN_SD = ',round(sd(TRAIN_RES),2)),
paste0('TEST_MEAN = ',round(mean(TEST_RES),2)),
paste0('TEST_SD = ',round(sd(TEST_RES),2))
)
)
p5 = ggplot()+
geom_density(data = df_res,aes(x = RES, y = ..count..,fill = Type))+
geom_text(data = label_df_res,aes(x = x, y = y, label = labels),
color = 'black',fontface= 'bold',size = 5,hjust= 1)+
theme_bw()+
theme(
axis.text = element_text(face = 'bold',color = 'black',size = 12),
axis.title = element_text(face = 'bold',color = 'black',size = 14, hjust = .5),
legend.text = element_text(face = 'bold',color = 'black',size = 12),
legend.title = element_blank(),
legend.position = 'bottom',
legend.direction = 'horizontal'
)+
xlab('RES')+
ylab('Count')
png('plot5.png',
height = 20,
width = 20,
units = 'cm',
res= 800)
print(p5)
dev.off()
图5 模型残差分析
10. 总结
本篇主要解决了以下问题：
1. 如何在R中构建核支持向量机模型？
2. 如何对其预测结果进行定量化评估？
11. 本篇所使用的软件包（没有的需要通过install.packages('包名')进行安装）
library(reshape2)
library(ggplot2)
library(e1071)
library(hydroGOF)
12. 致谢
首先，感谢大家的持续关注，小编会继续努力，持续更新下去的！
大家如果觉得有用，还麻烦大家转发点赞加关注哈，也可以扩散到朋友圈，多谢大家啦～
大家如果在使用本文代码的过程有遇到问题的，可以留言评论，也可以私信我哈~~
最后，祝大家即将到来的元旦节快乐！
小编联系方式
