工业产量分析与预测

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。

相关文章

Python的Darts库实现时间序列预测

Python的Darts库实现时间序列预测

《Python的Darts库实现时间序列预测》Darts一个集统计、机器学习与深度学习模型于一体的Python时间序列预测库,本文主要介绍了Python的Darts库实现时间序列预测,感兴趣的可以了解... 目录目录一、什么是 Darts？二、安装与基本配置安装 Darts导入基础模块三、时间序列数据结构与
阅读更多...
Android 缓存日志Logcat导出与分析最佳实践

Android 缓存日志Logcat导出与分析最佳实践

《Android缓存日志Logcat导出与分析最佳实践》本文全面介绍AndroidLogcat缓存日志的导出与分析方法,涵盖按进程、缓冲区类型及日志级别过滤,自动化工具使用,常见问题解决方案和最佳实... 目录android 缓存日志（Logcat）导出与分析全攻略为什么要导出缓存日志？按需过滤导出1. 按
阅读更多...
Linux中的HTTPS协议原理分析

Linux中的HTTPS协议原理分析

《Linux中的HTTPS协议原理分析》文章解释了HTTPS的必要性：HTTP明文传输易被篡改和劫持,HTTPS通过非对称加密协商对称密钥、CA证书认证和混合加密机制,有效防范中间人攻击,保障通信安全... 目录一、什么是加密和解密？二、为什么需要加密？三、常见的加密方式3.1 对称加密3.2非对称加密四、
阅读更多...
MySQL中读写分离方案对比分析与选型建议

MySQL中读写分离方案对比分析与选型建议

《MySQL中读写分离方案对比分析与选型建议》MySQL读写分离是提升数据库可用性和性能的常见手段,本文将围绕现实生产环境中常见的几种读写分离模式进行系统对比,希望对大家有所帮助... 目录一、问题背景介绍二、多种解决方案对比2.1 原生mysql主从复制2.2 Proxy层中间件：ProxySQL2.3
阅读更多...
python使用Akshare与Streamlit实现股票估值分析教程（图文代码）

python使用Akshare与Streamlit实现股票估值分析教程（图文代码）

《python使用Akshare与Streamlit实现股票估值分析教程（图文代码）》入职测试中的一道题,要求:从Akshare下载某一个股票近十年的财务报表包括,资产负债表,利润表,现金流量表,保存... 目录一、前言二、核心知识点梳理1、Akshare数据获取2、Pandas数据处理3、Matplotl
阅读更多...
python panda库从基础到高级操作分析

python panda库从基础到高级操作分析

《pythonpanda库从基础到高级操作分析》本文介绍了Pandas库的核心功能,包括处理结构化数据的Series和DataFrame数据结构,数据读取、清洗、分组聚合、合并、时间序列分析及大数据... 目录1. Pandas 概述2. 基本操作：数据读取与查看3. 索引操作：精准定位数据4. Group
阅读更多...
MySQL中EXISTS与IN用法使用与对比分析

MySQL中EXISTS与IN用法使用与对比分析

《MySQL中EXISTS与IN用法使用与对比分析》在MySQL中,EXISTS和IN都用于子查询中根据另一个查询的结果来过滤主查询的记录,本文将基于工作原理、效率和应用场景进行全面对比... 目录一、基本用法详解1. IN 运算符2. EXISTS 运算符二、EXISTS 与 IN 的选择策略三、性能对比
阅读更多...
MySQL 内存使用率常用分析语句

MySQL 内存使用率常用分析语句

《MySQL内存使用率常用分析语句》用户整理了MySQL内存占用过高的分析方法,涵盖操作系统层确认及数据库层bufferpool、内存模块差值、线程状态、performance_schema性能数据... 目录一、 OS层二、 DB层1. 全局情况2. 内存占js用详情最近连续遇到mysql内存占用过高导致
阅读更多...
深度解析Nginx日志分析与499状态码问题解决

深度解析Nginx日志分析与499状态码问题解决

《深度解析Nginx日志分析与499状态码问题解决》在Web服务器运维和性能优化过程中,Nginx日志是排查问题的重要依据,本文将围绕Nginx日志分析、499状态码的成因、排查方法及解决方案展开讨论... 目录前言1. Nginx日志基础1.1 Nginx日志存放位置1.2 Nginx日志格式2. 499
阅读更多...
Olingo分析和实践之EDM 辅助序列化器详解(最佳实践)

Olingo分析和实践之EDM 辅助序列化器详解(最佳实践)

《Olingo分析和实践之EDM辅助序列化器详解(最佳实践)》EDM辅助序列化器是ApacheOlingoOData框架中无需完整EDM模型的智能序列化工具,通过运行时类型推断实现灵活数据转换,适用... 目录概念与定义什么是 EDM 辅助序列化器？核心概念设计目标核心特点1. EDM 信息可选2. 智能类
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2