工业产量分析与预测

最新推荐文章于 2023-08-22 15:20:44 发布
最新推荐文章于 2023-08-22 15:20:44 发布
阅读量1.1k 收藏 11
点赞数 2
分类专栏： R语言数据分析与挖掘 文章标签： 数据分析 工业产量 预测 聚类 可视化
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_43180762/article/details/118388925
版权

工业产量分析与预测

导入数据与数据概况

# 安装库专用

# 通过如下命令设定镜像
options(repos = 'http://mirrors.ustc.edu.cn/CRAN/')
# 查看镜像是否修改
getOption('repos')
# 尝试下载R包
#若有需要，进行安装
#install.packages('plot3D')
#install.packages("reshape2")
#install.packages("pheatmap")

‘http://mirrors.ustc.edu.cn/CRAN/’

#设置工作路径
setwd("D:/LengPY")

#导入数据
library(readxl)
yhyear<-read_excel("yhdata.xlsx",sheet=1)
yhshengfen<-read_excel("yhdata.xlsx",sheet=2)

### 年份数据
head(yhyear)
A tibble: 6 × 4
年份钢材汽车天然气
<chr><dbl><dbl><dbl>
200537771.14 570.49493.20
200646893.36 727.89585.53
200756560.87 888.89692.40
200860460.29 930.59802.99
200969405.401379.53852.69
201080276.581826.53957.91
### 分省份数据
head(yhshengfen)
A tibble: 6 × 4
地区钢材汽车天然气
<chr><dbl><dbl><dbl>
北京 170.7077164.018629.73
天津 5454.9542104.149134.90
河北 28409.6316105.0836 5.84
山西 5594.2450 6.578064.62
内蒙古 2565.0281 2.897722.07
辽宁 7254.4315 79.1603 6.18

数据描述统计

summary(yhshengfen)

     地区                钢材            汽车             天然气       
 Length:31          Min.   :    0   Min.   :  0.000   Min.   :  0.000  
 Class :character   1st Qu.: 1037   1st Qu.:  6.135   1st Qu.:  0.185  
 Mode  :character   Median : 2565   Median : 61.858   Median :  5.840  
                    Mean   : 3886   Mean   : 82.828   Mean   : 56.829  
                    3rd Qu.: 3804   3rd Qu.:104.616   3rd Qu.: 40.280  
                    Max.   :28410   Max.   :311.972   Max.   :473.420  

summary(yhyear)

     年份                钢材             汽车            天然气      
 Length:15          Min.   : 37771   Min.   : 570.5   Min.   : 493.2  
 Class :character   1st Qu.: 64933   1st Qu.:1155.1   1st Qu.: 827.8  
 Mode  :character   Median : 95578   Median :1927.6   Median :1106.1  
                    Mean   : 86863   Mean   :1879.5   Mean   :1107.5  
                    3rd Qu.:106507   3rd Qu.:2509.0   3rd Qu.:1357.4  
                    Max.   :120457   Max.   :2901.8   Max.   :1761.7  

library(reshape2)
plotdata<-melt(yhyear,id.vars=c("年份"),variable.name="type",value.name="number")
plotdata$年份<-as.numeric(plotdata$年份)
head(plotdata)

Warning message:
"package 'reshape2' was built under R version 4.0.5"
A data.frame: 6 × 3
年份typenumber
<dbl><fct><dbl>
12005钢材37771.14
22006钢材46893.36
32007钢材56560.87
42008钢材60460.29
52009钢材69405.40
62010钢材80276.58

时序数据进行分析

#数据变形
library(reshape2)
mydata<-melt(yhyear,id.vars=c("年份"),variable.name="type",value.name="number")

head(mydata)
A data.frame: 6 × 3
年份typenumber
<chr><fct><dbl>
12005钢材37771.14
22006钢材46893.36
32007钢材56560.87
42008钢材60460.29
52009钢材69405.40
62010钢材80276.58
#绘制不同点形状的折线图
cardata<-mydata[which(mydata$type=='钢材'),]
library(ggplot2)
p <- ggplot(cardata,aes(x=年份,y=number,,shape=type,colour=type,group=type,fill=type)) +
			geom_line(size =0.8)+
            geom_point(mapping = NULL,data = NULL,size = 5,na.rm = FALSE,show.legend = NA,inherit.aes = TRUE)+
            geom_point(shape="\u2620")

p

Warning message:
"package 'ggplot2' was built under R version 4.0.4"

在这里插入图片描述

#绘制不同点形状的折线图
cardata<-mydata[which(mydata$type!='钢材'),]
library(ggplot2)
p <- ggplot(cardata,aes(x=年份,y=number,,shape=type,colour=type,group=type,fill=type)) +
			geom_line(size =0.8)+
            geom_point(mapping = NULL,data = NULL,size = 5,na.rm = FALSE,show.legend = NA,inherit.aes = TRUE)+
            geom_point(shape="\u2620")

p

在这里插入图片描述

尝试建立时间序列预测模型

钢材

library(ggplot2)
theme_set(theme_bw(base_family = "STKaiti"))
library(dplyr)
library(tidyr)
library(zoo)
library(tseries)
library(ggfortify)
library(gridExtra)
library(forecast)

​

# ## 转化为时间序列数据
gangcai_ts<-ts(yhyear$钢材)
## Ljung-Box 检验
Box.test(gangcai_ts,type ="Ljung-Box")

​

	Box-Ljung test

data:  gangcai_ts
X-squared = 10.972, df = 1, p-value = 0.000925

## 分析序列的自相关系数和偏自相关系数确定参数p和q
p1 <- autoplot(acf(gangcai_ts,lag.max = 30,plot = F))+
  ggtitle("序列自相关图")
p2 <- autoplot(pacf(gangcai_ts,lag.max = 30,plot = F))+
  ggtitle("序列偏自相关图")
gridExtra::grid.arrange(p1,p2,nrow=2)

在这里插入图片描述

auto.arima(gangcai_ts)

Series: gangcai_ts ARIMA(0,1,0) with drift Coefficients:         drift      5906.129s.e.  1441.853sigma^2 estimated as 31343862:  log likelihood=-140.17AIC=284.34   AICc=285.43   BIC=285.62

## 对数据建立ARIMA(1,1,0)模型，并预测后面的数据ARIMAmod1 <- arima(gangcai_ts,order = c(1,1,0))summary(ARIMAmod1)

Call:arima(x = gangcai_ts, order = c(1, 1, 0))Coefficients:         ar1      0.6761s.e.  0.1941sigma^2 estimated as 34933933:  log likelihood = -141.75,  aic = 287.51Training set error measures:                   ME     RMSE      MAE      MPE     MAPE      MASE       ACF1Training set 1951.221 5710.087 4604.412 2.718953 5.479596 0.6374903 -0.2117251

## 对拟合残差进行白噪声检验Box.test(ARIMAmod1$residuals,type ="Ljung-Box")

​

	Box-Ljung testdata:  ARIMAmod1$residualsX-squared = 0.8165, df = 1, p-value = 0.3662

## 可视化模型未来的预测值par(family = "STKaiti")plot(forecast(ARIMAmod1,h=5))

Warning message in title(main = main, xlab = xlab, ylab = ylab, ...):"Windows字体数据库里没有这样的字体系列"Warning message in axis(1, ...):"Windows字体数据库里没有这样的字体系列"Warning message in axis(2, ...):"Windows字体数据库里没有这样的字体系列"Warning message in axis(2, ...):"Windows字体数据库里没有这样的字体系列"Warning message in axis(2, ...):"Windows字体数据库里没有这样的字体系列"Warning message in axis(2, ...):"Windows字体数据库里没有这样的字体系列"Warning message in axis(2, ...):"Windows字体数据库里没有这样的字体系列"

在这里插入图片描述

# z这数据预测效果不太好，应该尝试多元方法

forecast(ARIMAmod1)

   Point Forecast     Lo 80    Hi 80     Lo 95    Hi 9516       125304.1 117729.53 132878.7 113719.78 136888.517       128581.2 113797.67 143364.7 105971.74 151190.718       130796.7 108896.36 152697.1  97303.01 164290.519       132294.6 103627.21 160962.0  88451.61 176137.620       133307.3  98306.37 168308.2  79778.01 186836.521       133991.9  93096.68 174887.2  71448.06 196535.822       134454.8  88077.02 180832.6  63526.11 205383.523       134767.7  83279.78 186255.7  56023.72 213511.724       134979.3  78711.94 191246.6  48925.81 221032.825       135122.3  74366.95 195877.7  42205.00 228039.6

先对汽车进行分析

# ## 转化为时间序列数据trq_ts<-ts(yhyear$汽车)## Ljung-Box 检验Box.test(trq_ts,type ="Ljung-Box")

​

	Box-Ljung testdata:  trq_tsX-squared = 12.75, df = 1, p-value = 0.0003561

p-value =0.001106<0.05 说明该序列为非随机数据

library(ggfortify)library(gridExtra)library(forecast)

autoplot(trq_ts)+ggtitle("序列变化趋势")

Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"

在这里插入图片描述

## 白噪声检验Box.test(trq_ts,type ="Ljung-Box")

​

	Box-Ljung testdata:  trq_tsX-squared = 12.75, df = 1, p-value = 0.0003561

p-value = 0.001106 ,说明不是白噪声

## 平稳性检验，单位根检验adf.test(trq_ts)

Warning message in adf.test(trq_ts):"p-value greater than printed p-value"




	Augmented Dickey-Fuller Testdata:  trq_tsDickey-Fuller = 0.33839, Lag order = 2, p-value = 0.99alternative hypothesis: stationary

p-value = 0.4931,说明数据是不平稳的

## 分析序列的自相关系数和偏自相关系数确定参数p和qp1 <- autoplot(acf(trq_ts,lag.max = 30,plot = F))+  ggtitle("序列自相关图")p2 <- autoplot(pacf(trq_ts,lag.max = 30,plot = F))+  ggtitle("序列偏自相关图")gridExtra::grid.arrange(p1,p2,nrow=2)

在这里插入图片描述

偏自相关图3阶后截尾，可以认为p的取值为1左右，

自相关图5阶后截尾，可以认为q的取值为0左右，???不确定

auto.arima(trq_ts)

Series: trq_ts ARIMA(0,1,0) with drift Coefficients:         drift      142.6560s.e.   49.9215sigma^2 estimated as 37574:  log likelihood=-93.08AIC=190.17   AICc=191.26   BIC=191.45

可知ARIMA(0,1,0)较好

## 对数据建立ARIMA(0,1,0)模型，并预测后面的数据ARIMAmod <- arima(trq_ts,order = c(0,1,0))summary(ARIMAmod)

Call:arima(x = trq_ts, order = c(0, 1, 0))

​

sigma^2 estimated as 55241:  log likelihood = -96.3,  aic = 194.6Training set error measures:                   ME     RMSE      MAE      MPE     MAPE      MASE      ACF1Training set 133.1836 227.0638 177.7351 8.847576 10.53488 0.9335331 0.2076724

由上表可知模型效果情况

## 对拟合残差进行白噪声检验Box.test(ARIMAmod$residuals,type ="Ljung-Box")

​

	Box-Ljung testdata:  ARIMAmod$residualsX-squared = 0.019715, df = 1, p-value = 0.8883

p-value = 0.3503 ,说明是白噪声

## 可视化模型未来的预测值par(family = "STKaiti")plot(forecast(ARIMAmod,h=5))

在这里插入图片描述

forecast(ARIMAmod)

   Point Forecast     Lo 80    Hi 80      Lo 95    Hi 9516       2448.930 2200.4523 2697.408 2068.91609 2828.94417       2383.369 1924.5864 2842.152 1681.72152 3085.01618       2347.171 1696.4904 2997.852 1352.04095 3342.30119       2327.185 1504.6088 3149.762 1069.16309 3585.20720       2316.151 1339.5995 3292.702  822.64448 3809.65721       2310.058 1194.5963 3425.520  604.10650 4016.01022       2306.694 1064.7215 3548.667  407.26076 4206.12823       2304.837  946.5300 3663.144  227.48573 4382.18824       2303.812  837.5700 3770.053   61.38863 4546.23525       2303.245  736.0695 3870.421  -93.54334 4700.034

试试灰色预测

都是用15年数据预测后五年，评估是看15年误差情况，后五个就是预测值

灰色预测汽车

yhyear
A tibble: 15 × 4
年份钢材汽车天然气
<chr><dbl><dbl><dbl>
2005 37771.14 570.490 493.200
2006 46893.36 727.890 585.530
2007 56560.87 888.890 692.400
2008 60460.29 930.590 802.990
2009 69405.401379.530 852.690
2010 80276.581826.530 957.910
2011 88619.571841.6401053.370
2012 95577.831927.6201106.080
2013108200.542212.0901208.580
2014112513.122372.5201301.570
2015103468.412450.3501346.097
2016104813.452811.9101368.650
2017104642.052901.8101480.350
2018113287.332782.7401601.590
2019120456.942567.6741761.740
#灰色预测x<-yhyear$汽车gm11<-function(x,k){#x为行向量数据#做一次累加n<-length(x)x1<-numeric(n);for(i in 1:n){   x1[i]<-sum(x[1:i]);}#x1的均值数列z1<-numeric(n)m<-n-1for(j in 1:m){   z1[j+1]<-(0.5*x1[j+1]+0.5*x1[j])}Yn=t(t(x[2:n]))B<-matrix(1,nrow=n-1,ncol=2)B[,1]<-t(t(-z1[2:n]))#solve(M)求M的逆#最小二乘法求解参数列u<-solve(t(B)%*%B)%*%t(B)%*%Yn;a<-u[1];b<-u[2];#预测x2<-numeric(k);x2[1]<-x[1];for(i in 1:k-1){   x2[1+i]=(x[1]-b/a)*exp(-a*i)+b/a;}x2=c(0,x2);#还原数据y=diff(x2);y}#调用函数gm11(x,20)
  1. 570.49
  2. 1151.7064043779
  3. 1245.64577762778
  4. 1347.2473517763
  5. 1457.13609716952
  6. 1575.98796009879
  7. 1704.53402067315
  8. 1843.56498983028
  9. 1993.93607314784
  10. 2156.57223137344
  11. 2332.47387003174
  12. 2522.72299310652
  13. 2728.48985865025
  14. 2951.04017726093
  15. 3191.74289770544
  16. 3452.07862758024
  17. 3733.64874080659
  18. 4038.18522798192
  19. 4367.56135017904
  20. 4723.80316172634
##qichecar<-data.frame(year=2005:2019,pre=gm11(x,15),yhyear$汽车)

head(car)
A data.frame: 6 × 3
yearpreyhyear.汽车
<int><dbl><dbl>
12005 570.490 570.49
220061151.706 727.89
320071245.646 888.89
420081347.247 930.59
520091457.1361379.53
620101575.9881826.53
library(Metrics)library(ROCR)#全部数据sprintf("MAPE误差为: %f",mape(car$pre,car$yhyear.汽车))

Warning message:"package 'Metrics' was built under R version 4.0.4"Attaching package: 'Metrics'

​

The following object is masked from 'package:forecast':    accuracy

​

Warning message:"package 'ROCR' was built under R version 4.0.4"

‘MAPE误差为: 0.132848’

plotdata <- data.frame(X =car$year,Y1 = car$yhyear.汽车,rfpre =car$pre)plotdata <- gather(plotdata,key="model",value="value",c(-X))## 可视化模型的预测误差ggplot(plotdata,aes(x = X,y = value))+  geom_line(aes(linetype = model,colour = model),size = 0.8)+  theme(legend.position = c(0.1,0.8),        plot.title = element_text(hjust = 0.5))+  ggtitle("预测效果可视化")

在这里插入图片描述

灰色预测天然气

#灰色预测x<-yhyear$天然气gm11<-function(x,k){#x为行向量数据#做一次累加n<-length(x)x1<-numeric(n);for(i in 1:n){   x1[i]<-sum(x[1:i]);}#x1的均值数列z1<-numeric(n)m<-n-1for(j in 1:m){   z1[j+1]<-(0.5*x1[j+1]+0.5*x1[j])}Yn=t(t(x[2:n]))B<-matrix(1,nrow=n-1,ncol=2)B[,1]<-t(t(-z1[2:n]))#solve(M)求M的逆#最小二乘法求解参数列u<-solve(t(B)%*%B)%*%t(B)%*%Yn;a<-u[1];b<-u[2];#预测x2<-numeric(k);x2[1]<-x[1];for(i in 1:k-1){   x2[1+i]=(x[1]-b/a)*exp(-a*i)+b/a;}x2=c(0,x2);#还原数据y=diff(x2);y}#调用函数gm11(x,20)
  1. 493.200000000001
  2. 698.574451446175
  3. 750.082742319664
  4. 805.388916186466
  5. 864.773003988903
  6. 928.535684311362
  7. 996.999805801801
  8. 1070.51202184655
  9. 1149.44454577538
  10. 1234.19703548378
  11. 1325.19861701499
  12. 1422.91005734758
  13. 1527.82609739019
  14. 1640.47795699596
  15. 1761.43602468017
  16. 1891.31274566023
  17. 2030.76572284046
  18. 2180.50104644332
  19. 2341.27686914573
  20. 2513.90724482249
##天然气tairanqi<-data.frame(year=2005:2019,pre=gm11(x,15),yhyear$天然气)

library(Metrics)library(ROCR)#全部数据sprintf("MAPE误差为: %f",mape(tairanqi$pre,tairanqi$yhyear.天然气))

‘MAPE误差为: 0.039464’

plotdata <- data.frame(X =tairanqi$year,Y1 = tairanqi$yhyear.天然气,rfpre =tairanqi$pre)plotdata <- gather(plotdata,key="model",value="value",c(-X))## 可视化模型的预测误差ggplot(plotdata,aes(x = X,y = value))+  geom_line(aes(linetype = model,colour = model),size = 0.8)+  theme(legend.position = c(0.1,0.8),        plot.title = element_text(hjust = 0.5))+  ggtitle("预测效果可视化")

Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"

在这里插入图片描述

灰色预测钢材

#灰色预测x<-yhyear$钢材gm11<-function(x,k){#x为行向量数据#做一次累加n<-length(x)x1<-numeric(n);for(i in 1:n){   x1[i]<-sum(x[1:i]);}#x1的均值数列z1<-numeric(n)m<-n-1for(j in 1:m){   z1[j+1]<-(0.5*x1[j+1]+0.5*x1[j])}Yn=t(t(x[2:n]))B<-matrix(1,nrow=n-1,ncol=2)B[,1]<-t(t(-z1[2:n]))#solve(M)求M的逆#最小二乘法求解参数列u<-solve(t(B)%*%B)%*%t(B)%*%Yn;a<-u[1];b<-u[2];#预测x2<-numeric(k);x2[1]<-x[1];for(i in 1:k-1){   x2[1+i]=(x[1]-b/a)*exp(-a*i)+b/a;}x2=c(0,x2);#还原数据y=diff(x2);y}#调用函数gm11(x,20)
  1. 37771.1399999999
  2. 62010.887920134
  3. 65519.7723170475
  4. 69227.2068415897
  5. 73144.4264473023
  6. 77283.3018172646
  7. 81656.3753368349
  8. 86276.8991019032
  9. 91158.8750778423
  10. 96317.0975308546
  11. 101767.197860289
  12. 107525.6919678
  13. 113610.030306885
  14. 120038.650764475
  15. 126831.034534831
  16. 134007.765155054
  17. 141590.590881134
  18. 149602.490593528
  19. 158067.743432017
  20. 167012.002370838
##天然气gangcai<-data.frame(year=2005:2019,pre=gm11(x,15),yhyear$钢材)

library(Metrics)library(ROCR)#全部数据sprintf("MAPE误差为: %f",mape(gangcai$pre,gangcai$yhyear.钢材))

‘MAPE误差为: 0.091504’

plotdata <- data.frame(X =gangcai$year,Y1 = gangcai$yhyear.钢材,rfpre =gangcai$pre)plotdata <- gather(plotdata,key="model",value="value",c(-X))## 可视化模型的预测误差ggplot(plotdata,aes(x = X,y = value))+  geom_line(aes(linetype = model,colour = model),size = 0.8)+  theme(legend.position = c(0.1,0.8),        plot.title = element_text(hjust = 0.5))+  ggtitle("预测效果可视化")

Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"

在这里插入图片描述

组合模型

结合时间序列和灰色预测，取个平均就行

探索相关关系

head(yhshengfen)
A tibble: 6 × 4
地区钢材汽车天然气
<chr><dbl><dbl><dbl>
北京 170.7077164.018629.73
天津 5454.9542104.149134.90
河北 28409.6316105.0836 5.84
山西 5594.2450 6.578064.62
内蒙古 2565.0281 2.897722.07
辽宁 7254.4315 79.1603 6.18
library(ggplot2)library(ggExtra)g <- ggplot(yhyear, aes(钢材, 汽车)) +geom_count() +geom_smooth(method="lm", se=F,color='red')ggMarginal(g, type = "histogram", fill="blue")ggMarginal(g, type = "boxplot", fill="transparent") ggMarginal(g, type = "density", fill="transparent")

在这里插入图片描述

表明两个产品有明显的线性关系

#(a) 二维散点与统计直方图library(ggplot2)library(RColorBrewer)library(ggpubr)ggscatterhist(  yhyear, x ='天然气', y = '汽车', shape=21,fill="#00AFBB",color = "black",size = 3, alpha = 1,  #palette = c("#00AFBB", "#E7B800", "#FC4E07"),  margin.params = list( fill="#00AFBB",color = "black", size = 0.2,alpha=1),  margin.plot =  "histogram",  legend = c(0.8,0.8),  ggtheme = theme_minimal())

在这里插入图片描述

同上，有线性关系

mcor<-cor(yhyear[,c(2:4)])library(corrplot) # 载入corrplot包corrplot(mcor)

在这里插入图片描述

mcor
A matrix: 3 × 3 of type dbl
钢材汽车天然气
钢材1.00000000.95439120.9575996
汽车0.95439121.00000000.9490545
天然气0.95759960.94905451.0000000

可发现，这三种产业都具有极高的相关性，可结合这几个行业的关系进行补充分析

各省份进行聚类分析及可视化

标准化后再聚类

cludata<-yhshengfen[,2:4]library(RSNNS)## 数据max-min归一化到0-1之间cludata[,1:3] <- normalizeData(cludata[,1:3],"0_1")

## 计算组内平方和  组间平方和tot_withinss <- vector()betweenss <- vector()for(ii in 1:8){  k1 <- kmeans(cludata[,c(1:3)],ii)  tot_withinss[ii] <- k1$tot.withinss  betweenss[ii] <- k1$betweenss}kmeanvalue <- data.frame(kk = 1:8,                         tot_withinss = tot_withinss,                         betweenss = betweenss)

library(ggplot2)library(gridExtra)library(ggdendro)library(cluster)library(ggfortify)p1 <- ggplot(kmeanvalue,aes(x = kk,y = tot_withinss))+  theme_bw()+  geom_point() + geom_line() +labs(y = "value") +  ggtitle("Total within-cluster sum of squares")+  theme(plot.title = element_text(hjust = 0.5))+  scale_x_continuous("kmean 聚类个数",kmeanvalue$kk)p2 <- ggplot(kmeanvalue,aes(x = kk,y = betweenss))+  theme_bw()+  geom_point() +geom_line() +labs(y = "value") +  ggtitle("The between-cluster sum of squares") +  theme(plot.title = element_text(hjust = 0.5))+  scale_x_continuous("kmean 聚类个数",kmeanvalue$kk)grid.arrange(p1,p2,nrow=2)

在这里插入图片描述

set.seed(245)k3 <- kmeans(cludata[,c(1:3)],5)summary(k3)

             Length Class  Mode   cluster      31     -none- numericcenters      15     -none- numerictotss         1     -none- numericwithinss      5     -none- numerictot.withinss  1     -none- numericbetweenss     1     -none- numericsize          5     -none- numericiter          1     -none- numericifault        1     -none- numeric

## 对聚类结果可视化clusplot(cludata[,c(1:3)],k3$cluster,main = "kmean cluster number=3")

在这里插入图片描述

## 可视化轮廓图，表示聚类效果sis1 <- silhouette(k3$cluster,dist(cludata[,c(1:3)],method = "euclidean"))plot(sis1,main = " kmean silhouette",     col = c("red", "green", "blue"))

在这里插入图片描述

聚类效果良好，有效地将各省份区分开，可结合标签和汇总表分析各类比的区别

#将标签写入cludata$bzcluster<-k3$clusterhead(cludata)
A tibble: 6 × 4
钢材汽车天然气bzcluster
<dbl><dbl><dbl><int>
0.0060087970.5257484980.062798362
0.1920107330.3338416060.073718903
1.0000000000.3368370710.012335771
0.1969136770.0210852530.136496135
0.0902872700.0092883450.046618225
0.2553511280.2537420070.013053953
cludata$地区<-yhshengfen$地区cludata
A tibble: 31 × 5
钢材汽车天然气bzcluster地区
<dbl><dbl><dbl><int><chr>
0.0060087970.52574849766.279836e-022北京
0.1920107330.33384160617.371890e-023天津
1.0000000000.33683707111.233577e-021河北
0.1969136770.02108525261.364961e-015山西
0.0902872700.00928834554.661822e-025内蒙古
0.2553511280.25374200731.305395e-023辽宁
0.0543562400.92676288482.181995e-022吉林
0.0275254940.06056448739.644713e-025黑龙江
0.0640517770.88118181272.634025e-022上海
0.5002320510.26867413572.574881e-023江苏
0.1220800760.31795971172.112289e-053浙江
0.1111724530.24880437834.478053e-033安徽
0.1315630050.05179349660.000000e+005福建 
aggregate(cludata[,c(1:3)],list(cludata$bzcluster),mean)
A data.frame: 5 × 4
Group.1钢材汽车天然气
<int><dbl><dbl><dbl>
11.000000000.336837070.01233577
20.088943380.773046330.05976018
30.186694230.266826790.02716193
40.078772250.129773890.88490840
50.054825850.017978880.03854544
yhshengfen$bzcluster<-k3$cluster

head(yhshengfen)
A tibble: 6 × 5
地区钢材汽车天然气bzcluster
<chr><dbl><dbl><dbl><int>
北京 170.7077164.018629.732
天津 5454.9542104.149134.903
河北 28409.6316105.0836 5.841
山西 5594.2450 6.578064.625
内蒙古 2565.0281 2.897722.075
辽宁 7254.4315 79.1603 6.183

不同类别间，各指标平均值，可得到各类别侧重产业，结合现实分析

#各类比分类统计aggregate(yhshengfen[,c(2:4)],list(yhshengfen$bzcluster),mean)
A data.frame: 5 × 4
Group.1钢材汽车天然气
<int><dbl><dbl><dbl>
128409.632105.08360 5.84000
2 2526.849241.16850 28.29167
3 5303.914 83.24238 12.85900
4 2237.891 40.48577418.93333
5 1557.582 5.60890 18.24818
head(yhshengfen)
A tibble: 6 × 5
地区钢材汽车天然气bzcluster
<chr><dbl><dbl><dbl><int>
北京 170.7077164.018629.732
天津 5454.9542104.149134.903
河北 28409.6316105.0836 5.841
山西 5594.2450 6.578064.625
内蒙古 2565.0281 2.897722.075
辽宁 7254.4315 79.1603 6.183

不同类别省份间，各产品产量分布情况。

library(treemap)treemap(yhshengfen,index = c("bzcluster","地区"),vSize = "钢材",        vColor = "汽车",type="value",palette="RdYlGn",        title = "不同分类下各地区钢材与汽车产量关系",fontfamily.title = "STKaiti",        title.legend = "汽车产量(万辆)",fontfamily.legend="STKaiti")

Warning message:"package 'treemap' was built under R version 4.0.4"Warning message in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"

在这里插入图片描述

library(treemap)treemap(yhshengfen,index = c("bzcluster","地区"),vSize = "天然气",        vColor = "汽车",type="value",palette="RdYlGn",        title = "不同分类下各地区天然气与汽车产量",fontfamily.title = "STKaiti",        title.legend = "汽车产量(万辆)",fontfamily.legend="STKaiti")

Warning message in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"

在这里插入图片描述

 ## 可视化各地区产量分布 ggplot(yhshengfen,aes(x=reorder(地区,汽车),y=汽车))+  theme_bw(base_family = "STKaiti")+  geom_bar(aes(fill=汽车),stat = "identity",show.legend = F)+  coord_flip()+  scale_fill_gradient(low = "#56B1F7", high = "#00C3C2")+  labs(x="地区",y="汽车产量",title="不同地区汽车产量")+  theme(axis.text.x = element_text(vjust = 0.5),        plot.title = element_text(hjust = 0.5))

Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :"Windows字体数据库里没有这样的字体系列"

在这里插入图片描述

俱往矣`
关注
  • 2
    点赞
  • 踩
  • 11
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
专栏目录
R ggplot Rmarkdown 输出乱码
Cynthia Pink Generation
11-22 6060
大面积出现如下  ## Warning in grid.Call(L_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : ## 'mbcsToSbcs'里转换'11鏈 21 06:00'出错：代替了dot ## Warning in grid.Call.graphics(L_text, as.graphicsAnnot(x$l
随笔7：R语言绘图黑白印刷风格芯片包
ForvaluesMyself的博客
01-18 909
R语言绘图黑白印刷风格芯片包 1.基本说明 1.1本包的目标 让使用者将注意力集中在参数设定与调整上，而不是代码使用的细节，类似于一个芯片封装的过程。对几个常用的绘图函数进行了封装，并调整为印刷风格供使用者做课题或写论文参考。 1.2功能特色 设定好参数后，只需要一行代码便能绘制出印刷风格的R图像。 1.3数据说明 本包内置数据集均为随机或手动生成，输入以下命令能够看到对应绘图函数的详细说明。 help(acchist) #堆积柱状图黑白印刷风格绘制 help(bonbon) ..
基于支持向量机的纺纱质量预测模型研究
01-15
<p>纱线的生产是一个多环节的复杂工业过程,其质量控制大多需要依赖领域专家的个人经验,为此提出一种基于支持向量机的纱线质量预测模型.探讨了模型选择与验证问题,并利用网格搜索法对模型参数进行了优化.试验结果表明,在小样本和噪音数据环境下,支持向量机模型仍能保持一定的预测精度,与人工神经网络模型相比,更适应于真实纺纱生产过程.</p>
工业生产流程的数据分析预测
08-29
工业故障智能预测：使用随机深林进行预测。含有：testA_suzhou，train_suzhou，测试A_20180117，测试B_20180117，训练_20180117等数据和模型代码
Python实现基于RF随机森林的AQI预测模型（完整源码和数据）
08-12
Python实现基于RF随机森林的AQI预测模型（完整源码和数据） Python实现基于RF随机森林的AQI预测模型（完整源码和数据） Python实现基于RF随机森林的AQI预测模型（完整源码和数据） 空气质量（air quality）的好坏反映了空气污染程度，它是依据空气中污染物浓度的高低来判断的。空气污染是一个复杂的现象，在特定时间和地点空气污染物浓度受到许多因素影响。来自固定和流动污染源的人为污染物排放大小是影响空气质量的最主要因素之一，其中包括车辆、船舶、飞机的尾气、工业企业生产排放、居民生活和取暖、垃圾焚烧等。城市的发展密度、地形地貌和气象等也是影响空气质量的重要因素。
Python实现基于MLR多元线性回归的AQI预测模型（完整源码和数据）
08-12
Python实现基于MLR多元线性回归的AQI预测模型（完整源码和数据） Python实现基于MLR多元线性回归的AQI预测模型（完整源码和数据） Python实现基于MLR多元线性回归的AQI预测模型（完整源码和数据） 空气质量（air quality）的好坏反映了空气污染程度，它是依据空气中污染物浓度的高低来判断的。空气污染是一个复杂的现象，在特定时间和地点空气污染物浓度受到许多因素影响。来自固定和流动污染源的人为污染物排放大小是影响空气质量的最主要因素之一，其中包括车辆、船舶、飞机的尾气、工业企业生产排放、居民生活和取暖、垃圾焚烧等。城市的发展密度、地形地貌和气象等也是影响空气质量的重要因素。
ERP沙盘模拟产能计算模型
09-18
ERP沙盘模拟产能计算模型，辅助计算产能，帮助用户选择销售订单
Error in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : 字体类别出错In addition: T
weixin_54434521的博客
08-22 2553
font_import("C:/Windows/Fonts")#选择yes。#再使用ggsave(）就不会报错。#导入windows中的字体库。#和Times相关的有下面这些。
【无标题】rstudio绘制的图形导出pdf不显示文字内容，导出png正常显示
weixin_47151413的博客
07-14 1490
PDF打不开处理办法
rstudio绘制的图形导出pdf不显示文字内容，导出png正常显示
yangqijia1的博客
10-10 4713
问题： 之前绘制一个热图的时候，发现导出pdf的时候居然不显示文字！但是导出png的时候又是是正常显示的。我就纳闷了，如图所示： csdn自带的编辑器不太会调整，各位见谅将就着看哈。 原因： （手动狗头）我发现我代码里面设置了字体，然鹅我系统里面没有这种字体。 pheatmap(fontfamily="Arial",fontsize = 6,fontface="italic") ### 这个不是完整的代码，是出问题的原因 Warning messages: 1: In grid.Call(C_text
Python实现基于GradientBoosting的AQI预测模型（完整源码和数据）
08-12
Python实现基于GradientBoosting的AQI预测模型（完整源码和数据） Python实现基于GradientBoosting的AQI预测模型（完整源码和数据） Python实现基于GradientBoosting的AQI预测模型（完整源码和数据） 空气质量（air quality）的好坏反映了空气污染程度，它是依据空气中污染物浓度的高低来判断的。空气污染是一个复杂的现象，在特定时间和地点空气污染物浓度受到许多因素影响。来自固定和流动污染源的人为污染物排放大小是影响空气质量的最主要因素之一，其中包括车辆、船舶、飞机的尾气、工业企业生产排放、居民生活和取暖、垃圾焚烧等。城市的发展密度、地形地貌和气象等也是影响空气质量的重要因素。
论文研究-能源工业生产量及主要能源消费量的分析与预测.pdf
09-20
论文研究-能源工业生产量及主要能源消费量的分析与预测.pdf, 本文利用文[1]中提出的关联度识别方法, 以云南省为例, 分析云南省主要能源(煤炭、石油、天然气、电力)消费...
2024-2030中国工业碳刷市场现状研究分析与发展前景预测报告.docx
最新发布
05-13
【重点】本报告研究全球与中国的产能、产量、销量、销售额、价格及未来趋势。重点分析全球与中国市场的主要厂商产品特点、产品规格、价格、销量、销售收入及全球和中国市场主要生产商的市场份额。 此外针对行业产品...
2022-2028全球与中国工业温度控制器市场现状及未来发展趋势.doc
01-10
本报告研究全球与中国市场工业温度控制器的产能、产量、销量、销售额、价格及未来趋势。重点分析全球与中国市场的主要厂商产品特点、产品规格、价格、销量、销售收入及全球和中国市场主要生产商的市场份额。历史数据...
电子玻璃行业深度分析：预测2024年我国产量将增至6.5亿平方米（简报）.pdf
03-18
横跨如半导体产业链（半导体设备及零部件、半导体材料、集成电路、制造、封测、分立器件、传感器、光电器件）、光伏产业链（设备、硅料/硅片、...化、工业机器人、激光、工控、无人机）、食品药品、医疗器械、农业等...
R语言实现灰色预测模型gm11
热门推荐
yuanzhoulvpi的博客
09-05 1万+
不多说，直接上R语言代码 以后有时间再继续更新 第一个先计算关联度和关联系数 #定义成函数------------------------------------------- hsyc <- function(y, ro) { #这个函数是计算关联度和关联度系数 #初始化 (x <- y / y[ ,1]) #求绝对差序列 x0 <- x[1...
报错记录：Error in grid.Call(C_stringMetric, as.graphicsAnnot(x$label))
生信练习生
11-29 696
Error in grid.Call(C_stringMetric, as.graphicsAnnot(x$label))
ggplot2 保存图片字体错误问题
Gossie的博客
09-12 8043
用ggplot2作完图，输出图片或者用ggsave()保存图片时可能会碰到字体错误问题，如设置字体为Arial后： Error in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : 字体类别出错 一种解决方法是用 showtext 包输入需要的字体： library(showtext) font_add('Arial','/Library/Fonts/Arial.ttf') #加载字体，MAC 中字体库在 /Lib
工业蒸汽预测 RStudio
06-08
您需要对工业蒸汽的产生进行预测，可以使用R语言及其集成开发环境RStudio来完成。具体步骤如下： 1. 收集数据：收集有关工业蒸汽产生的数据，包括温度、湿度、压力等因素的变化情况以及相应的蒸汽产量。 2. 数据清洗：对数据进行清洗，包括缺失值、异常值的处理。 3. 数据可视化：使用R语言中的ggplot2等数据可视化工具对数据进行探索性分析，以便更好地理解数据特征和规律。 4. 建立模型：根据数据特征和规律，选择合适的机器学习算法建立预测模型。 5. 模型评估：对模型进行评估，包括准确度、召回率、F1值等指标。 6. 模型部署：将预测模型部署到生产环境中，实现实时预测。 以上是一个大致的流程，具体细节需要根据实际情况进行调整。希望对您有所帮助。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值