工业产量分析与预测

2024-03-27 04:59
文章标签 分析 工业 产量 预测

本文主要是介绍工业产量分析与预测，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

工业产量分析与预测

    • 导入数据与数据概况
    • 数据描述统计
    • 时序数据进行分析
    • 尝试建立时间序列预测模型
    • 先对汽车进行分析
    • 试试灰色预测
      • 灰色预测汽车
      • 灰色预测天然气
      • 灰色预测钢材
    • 组合模型
      • 探索相关关系
    • 各省份进行聚类分析及可视化
      • 标准化后再聚类

导入数据与数据概况

# 安装库专用# 通过如下命令设定镜像
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 testdata:  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字体数据库里没有这样的字体系列"

在这里插入图片描述

这篇关于工业产量分析与预测的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/851008。

相关文章

SpringBoot中六种批量更新Mysql的方式效率对比分析

SpringBoot中六种批量更新Mysql的方式效率对比分析

《SpringBoot中六种批量更新Mysql的方式效率对比分析》文章比较了MySQL大数据量批量更新的多种方法,指出REPLACEINTO和ONDUPLICATEKEY效率最高但存在数据风险,MyB... 目录效率比较测试结构数据库初始化测试数据批量修改方案第一种 for第二种 case when第三种
阅读更多...
解决1093 - You can‘t specify target table报错问题及原因分析

解决1093 - You can‘t specify target table报错问题及原因分析

《解决1093-Youcan‘tspecifytargettable报错问题及原因分析》MySQL1093错误因UPDATE/DELETE语句的FROM子句直接引用目标表或嵌套子查询导致,... 目录报js错原因分析具体原因解决办法方法一：使用临时表方法二：使用JOIN方法三：使用EXISTS示例总结报错原
阅读更多...
MySQL中的LENGTH()函数用法详解与实例分析

MySQL中的LENGTH()函数用法详解与实例分析

《MySQL中的LENGTH()函数用法详解与实例分析》MySQLLENGTH()函数用于计算字符串的字节长度,区别于CHAR_LENGTH()的字符长度,适用于多字节字符集（如UTF-8）的数据验证... 目录1. LENGTH()函数的基本语法2. LENGTH()函数的返回值2.1 示例1：计算字符串
阅读更多...
Android kotlin中 Channel 和 Flow 的区别和选择使用场景分析

Android kotlin中 Channel 和 Flow 的区别和选择使用场景分析

《Androidkotlin中Channel和Flow的区别和选择使用场景分析》Kotlin协程中,Flow是冷数据流,按需触发,适合响应式数据处理；Channel是热数据流,持续发送,支持... 目录一、基本概念界定FlowChannel二、核心特性对比数据生产触发条件生产与消费的关系背压处理机制生命周期
阅读更多...
怎样通过分析GC日志来定位Java进程的内存问题

怎样通过分析GC日志来定位Java进程的内存问题

《怎样通过分析GC日志来定位Java进程的内存问题》：本文主要介绍怎样通过分析GC日志来定位Java进程的内存问题,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录一、GC 日志基础配置1. 启用详细 GC 日志2. 不同收集器的日志格式二、关键指标与分析维度1.
阅读更多...
MySQL中的表连接原理分析

MySQL中的表连接原理分析

《MySQL中的表连接原理分析》：本文主要介绍MySQL中的表连接原理分析,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录1、背景2、环境3、表连接原理【1】驱动表和被驱动表【2】内连接【3】外连接【4编程】嵌套循环连接【5】join buffer4、总结1、背景
阅读更多...
python中Hash使用场景分析

python中Hash使用场景分析

《python中Hash使用场景分析》Python的hash()函数用于获取对象哈希值,常用于字典和集合,不可变类型可哈希,可变类型不可,常见算法包括除法、乘法、平方取中和随机数哈希,各有优缺点,需根... 目录python中的 Hash除法哈希算法乘法哈希算法平方取中法随机数哈希算法小结在Python中，
阅读更多...
Java Stream的distinct去重原理分析

Java Stream的distinct去重原理分析

《JavaStream的distinct去重原理分析》Javastream中的distinct方法用于去除流中的重复元素,它返回一个包含过滤后唯一元素的新流,该方法会根据元素的hashcode和eq... 目录一、distinct 的基础用法与核心特性二、distinct 的底层实现原理1. 顺序流中的去重
阅读更多...
关于MyISAM和InnoDB对比分析

关于MyISAM和InnoDB对比分析

《关于MyISAM和InnoDB对比分析》：本文主要介绍关于MyISAM和InnoDB对比分析,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录开篇：从交通规则看存储引擎选择理解存储引擎的基本概念技术原理对比1. 事务支持：ACID的守护者2. 锁机制：并发控制的艺
阅读更多...
MyBatis Plus 中 update_time 字段自动填充失效的原因分析及解决方案(最新整理)

MyBatis Plus 中 update_time 字段自动填充失效的原因分析及解决方案(最新整理)

《MyBatisPlus中update_time字段自动填充失效的原因分析及解决方案(最新整理)》在使用MyBatisPlus时,通常我们会在数据库表中设置create_time和update... 目录前言一、问题现象二、原因分析三、总结：常见原因与解决方法对照表四、推荐写法前言在使用 MyBATis
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2