某地区共享电动汽车交易量使用情况
摘要
随着人口的快速增长和经济科技的高速发展，汽车的发展水平也是水涨船高，共享电动汽车最近几年也开始变得流行起来。共享电动汽车的发展相比于共享单车，起步还是比较晚的，这主要一个是技术上的问题，一个是经济上的问题，只有这两方面都成熟了或者满足一定条件，这种“流行”才能转变真正的社会现象。
本文将以中国北方某城市的某年每日的共享电动汽车交易数据为例，结合当地的天气气候环境，再结合节假日等信息，分析这些因素是否是影响电动汽车的交易量的原因，这里将通过可视化图表对信息陈列，还要经过相关检验知识更进一步确定；接着我们将对于每日交易量数据进行一些检验，判定此时序数据是不是前后有相关性的，我们可以结合时间序列知识进行相应检验。
由于数据是时间序列的特点，所以我们还将对这些数据做未来的预测分析，我们将建立两种模型进行预测。一种是基于深度学习的LSTM模型，一种是基于机器学习的SVM模型。此两种模型底层原理都是一种回归类型的模型，独立建立好这两种模型后，我们再结合上面所提到的“是否节假日”、“天气气候”等因子综合考虑，作为模型优化的一个方面。由于此分析报告重点是基于实际数据做出一些有价值的分析，中间涉及的模型和方法重要理解、运用和优化，底层理论公式推导将不再涉及。
关键字：时间序列分析、描述性统计、因子分析、非参数检验、LSTM、SVM
目录
一：引言
1：背景
随着共享单车在全国各大城市迅速铺开，“共享经济”的概念迅速普及，共享汽车也随之悄然进入了人们的视野。2015年，被外界定义为中国共享汽车的元年。从2015年开始区域型试点共享汽车项目开始规模化扩张，从华东地区零星几家到大江南北涌现出200多家企业。
这些共享汽车平台也像共享单车的发展模式，率先在北京、上海和广州等大型城市布局，虽然各家平台投放车辆以及网点的数量有多有少，但均已在市场上引起一定反响。共享汽车的出现为人们的生活带来方便快捷，但新事物的诞生还是需要和社会实际情况不断的磨合相适应。为了想让共享电动汽车未来更加科学有效的发展，我们就要结合科学有效的分析方法为其保驾护航。
2：分析意义
由于开始各个公司的运营模式或者成本问题，在现实生活中存在很多种棘手问题，比如刚开始汽车到底投入多少呢？这一方面是为了公司节约成本，利益最大化，另一方面也不会给社会造成空间上的浪费，所以这投入多少合适就涉及很多因素。比如在实际中，节假日的到来，一家人出行很大可能都会乘坐私家车，所以周末和节假日就是个考虑因素。如果运用好分析好这些因素，那么很多公司将避免掉很多问题，比如有可能开始投入汽车过大，成本开销过大，导致最后押金返还不了的情况。所以，我们基于科学的研究方法，来探讨当前交易量的影响因素和预测未来交易量的走势，是非常重要的，是有很大商业和实际的价值。
二：数据清洗
我们知道，在不改变原始数据的绝大多数信息的情况下，更加“干净无噪声”的数据将对分析起到事倍功半的效果，所以数据清洗是非常重要的建模分析步骤！
1：数据类型的改变
图1
由上图1我们发现，此数据的一些值的属性并不属于标准的数值格式，这样的是不能进行数值运算的，为了充分利用数据信息，我们将做数据类型的转换。
1）：“最低最高温度”属性的处理
我们把两个极值温度取平均作为均温度来参加运算。
2）：“天气”属性的处理
“天气”列涉及的信息很多，但从实际角度出发，我们主要关心“下雨”与“非下雨”，这两种情况对我们日常行为有比较大的影响，无论是大雨还是小雨，都设为下雨天，其他的天气设为非下雨天，下雨天我们设值为1，非下雨天为0。
修改后的数据，如下展示：
图2
2：缺失值的处理
我们编写相关程序查找这些属性值是否存在缺失情况。
图3
从上图3发现，只有“订单量”数据存在缺失情况，缺失占总比不是非常高。
对于缺失值的处理，我们主要用数据的中位数填充，填充结果如下：
图4
3：异常值的处理
处理异常值也是个非常重要的数据清洗过程，如果数据存在明显的噪声，那么噪声对数据的影响可以说是“牵一发而动全身”的效果，在建模过程中，模型的效果会因为噪声的存在，改变训练方向，所以我们非常有必要去除一些明显噪声。
我们将采用综合法对数据进行异常值处理，即用数据拟合法再结合统计法对数据整体去除异常值，这种方法的优点一是效率比较高，二是能充分利用所有数据信息对数据做异常值处理。
由于自身数据的特点，在粗略了观察“周日_2周六节假日_1班_0”和“非雨_0雨1”两属性队形的数据没有异常分类值后，我们将只对“订单量”和“均温度”涉及的数据进行数据异常值处理。首先我们先直观感受一下这两组数据的分布，然后我们再给出异常值处理后的数据分布图，综合判定我们的异常值处理是合理。
图5
图6
从图5和图6我们可以大致看出，“订单量”数据缺失存在明显的异常值和缺失值，而“均温度”数据看上去是比较正常的，随着季节的变化有很明显的周期性。接下来，我们看看处理后的效果。
图8
图7
上图7和图8就是通过“拟合-统计法”得到的异常值处理效果，我们发现“订单量”的异常数据被明显的替换了，“均温度”数据整体上没有明显的异常值出现，替换数据的逻辑为：
当异常值大于拟合上限值时，我们有替换值：
如果
那么异常值
否则，这样做的目的是此时的异常值替换不超过拟合上限值，防止拟合上限值过于小，同理当异常值小于拟合下限值时，我们有替换值：
如果
那么异常值
否则，这样做的目的是此时的异常值替换不小于拟合下限值，防止拟合下限值过于大。
图8
上图8就是展现我们异常值处理的计算过程，至此，我们的数据清洗阶段算是彻底完成，接下来就是关于更细致的分析和建模，但更加“干净无噪声”的数据将对后面分析起到事倍功半的效果！
三：数据探索式分析
在文章开头，我们提到过，共享电动车的交易量也许跟气候和节假日有关系，所以，带着这种疑问，我们将通过统计分析和可视化图表，以及相关统计检验来证实我们的假设。
1：因子可视化分析
首先我们基于处理好的数据，把“周日”、“周六和法定节假日”、“工作日”等信息对应的所有订单交易量数据进行求平均，然后再以柱状图的形式呈现在如下图9中。
图9
图10
图11
同理，我们对“是否是雨天”和“均温度”做同样的处理，得到图10和图11所示的柱状图信息。从柱状图的信息，我们大致能观察出，这些因素对订单量的多少是有些影响的，最后我们再观察下订单交易量的具体区间密度分布情况，如下图展示。
图12
从上图12我们发现，一年的订单交易量，发生在每天200-400单的概率是比较大的，最起码都有100单以上，总体概率密度分布不是非常符合标准的正态分布。从描述性表格数据来看，更加证实了上述统计图给我们的直观感觉。
2：相关检验分析
由于上述得出的结论很大都是出于直观感觉，所以我们有必要进行相关统计检验，来证实这些因子的不同确实对订单量是有影响的。从实际订单交易量数据分布图来看，我们的数据确实不是太符合正态分布，所以我们首先考虑利用非参数检验方法来进行相关检验。
2.1：Kruskal-Wallis（简称克氏）检验
作为对样本分布没有太大要求的Kruskal-Wallis（简称克氏）检验，它是一个将两个独立样本Wilcoxon(Mann-Whitney)推广到3个或者更多组的检验，可以检验多组独立数据均值（或者分布）之间的差异性。接下来，我们通过编程实现上面的检验，检验结果如下：
图13
从图13的p值结果我们发现，各个因素的不同对订单的影响是非常显著的，尤其“均温度”这一影响，更加的显著。
由于订单量数据是基于时间的，所以可以称为是时间序列数据，因此我们有必要对此数据进行相关操作，比如白噪声检验、数据结构分解等等。
2.2：自相关性检验
白噪声检验又称为自相关性检验，由于是时序数据，很可能存在一个前后数据有关系的情况，一旦数据前后有关系，在做未来预测的时候，很有可能后面的预测值是前面的预测值的一个“右偏”结果，一旦发生“右偏”现象，其实预测是无效的，因为未来预测的数据已经在过去发生了，预测值滞后于过去值，没任何实际意义，所以我们会对此问题做出优化。
图14
上图14就是对数据进行滞后15阶的一个自相关性检验，由p值发现我们的数据并不是白噪声数据，是有前后关系的，可以进行时间序列建模，但是要防止“右偏”现象的发生。最后，我们再对此数据进行一个结构分解，当我们确定好数据分解周期后（周期设为14天），得到如下分解图，发现我们的数据集主要还是存在一个周期季节性，没有明显的趋势性。
图15
四：数据深度分析
这一章节主要是基于一些模型对时序数据做预测，我们将分别利用机器学习模型和深度学习模型进行预测分析，然后再结合因子分析的内容，作为我们优化模型的选择。
1：概述
机器学习(ML, Machine Learning)是一门多学科交叉专业，涵盖概率论知识，统计学知识，复杂算法知识，使用计算机模拟人类学习方式，并将现有内容进行知识结构划分来有效提高学习效率。另外机器学习包括深度学习。
支持向量机（Support Vector Machine, SVM）是一类按监督学习方式，其决策边界是对学习样本求解的最大边距超平面，SVM使用铰链损失函数计算经验风险并在求解系统中加入了正则化项以优化结构风险，是一个具有稀疏性和稳健性的分类器。SVM可以通过核方法进行非线性分类。
深度学习(DL, Deep Learning)是机器学习领域中一个新的研究方向，它被引入机器学习使其更接近于最初的目标——人工智能。深度学习是学习样本数据的内在规律和表示层次，这些学习过程中获得的信息对诸如文字，图像和声音等数据的解释有很大的帮助。
长短期记忆网络（LSTM，Long Short-Term Memory）是一种时间循环神经网络，是为了解决一般的RNN（循环神经网络）存在的长期依赖问题而专门设计出来的，所有的RNN都具有一种重复神经网络模块的链式形式。LSTM适合于处理和预测时间序列中间隔和延迟非常长的重要事件。
2：评价指标
由于我们的任务是建立回归模型来预测，相比于分类器的一些评价指标如准确率、ROC、召回率等等，此任务的评价指标只选必要的两种:
2.1:拟合优度系数
拟合优度（Goodness of Fit）是指回归直线对观测值的拟合程度。度量拟合优度的统计量是可决系数（亦称确定系数）R²。修正后的R²最大值为1，最小值为0。R²的值越接近1，说明回归直线对观测值的拟合程度越好；反之，R²的值越小，说明回归直线对观测值的拟合程度越差。未修正的R²值可以为负数，即拟合效果非常差！
2.2：均方误差
均方误差（mean-square error, MSE）是反映估计量与被估计量之间差异程度的一种度量。一般地，在样本量一定时，评价一个点估计的好坏标准使用的指标总是点估计与参数真值的距离的函数，最常用的函数是距离的平方，由于估计量具有随机性，可以对该函数求期望，这就是下式给出的均方误差：
自然，我们希望估计的均方误差越小越好。
3：利用SVM建模
SVM模型是经典的机器学习模型，我们将利用此模型建立合适的回归模型来对“订单量”时序数据做预测分析。
3.1：优化前的模型建立
模型未优化之前，我们仅利用时序本身的数据进行建模。首先我们要对数据集进行划分，由于是时序数据，考虑到样本数据不是非常大的客观事实，所以我们选择数据前80%的作为训练数据，后20%作为预测数据，然后我们还要确定划分时序数据的步长，根据数据量长度出发，我们选择步长为5开始划分数据集（当然这个参数是可以变的，也是优化模型的一个方向）， 再接着我们利用SVM模型框架，主要对四个参数进行选择，一个是合适的核函数，一个是合适的惩罚系数，一个是核函数阶数，最后一个是gamma值。
确定好上面的参数后，我们开始输入数据并做出相关的预测效果图，如下图16所示：
图16
从上图观察，我们发现无论是在训练集预测阶段，还是测试集预测阶段，预测数据都没有发生明显的“右偏”现象，说明我们的模型建立是有效的，接着我们呈现如下评价指标：
图17
上图17就是本章节我们开始提到的两种评价模型指标，我们发现尽管总体R²值不错，但是测试集R²值略小，一方面是由于数据集不多的原因，模型需要更多的数据来训练，另一方面我们也需要对模型进行深度优化。其中的MSE值，我们要通过不同建模结果作对比来评判。
3.2：优化后的模型建立
由开始文章所说，订单量是受一些因素影响的，比如天气和温度等等，那何必不把这些因素一起放到模型中，综合建模？所以这是个非常重要的优化思路，那么加了因子分析就意味着模型变得更加复杂，就意味着建模难度加大，但为了保证模型是因为加了因子才变得更加优越，我们必须保持一些参数不变，比如划分数据集步长，比如核函数的选择等等，完成了这些准备工作，我们将编写程序开始建模，并把建模预测效果图呈现如下：
图18
从上图18观察，我们发现无论是在训练集预测阶段，还是测试集预测阶段，预测数据完全没有“右偏”现象，说明我们的模型建立更是有效的，接着我们呈现如下评价指标：
图19
从上图19的评价指标我们发现，无论是R²值还是MSE值，都比优化前的数值提升不少，这说明我们优化的方向是非常科学非常有效的，但这是用模型的复杂度提升和运行程序耗时增加换来的。
4：利用LSTM网络建模
4.1：优化前的模型建立
我们先对数据进行“mean-std”标准化处理，然后确定好输入神经元个数，隐藏层层数，隐藏层神经元个数和输出神经元个数，激活函数和dropout比例，接着再设定预测步长，跟SVM建模一样也是选5，损失函数选择MSE，优化器为Adam，最后再设定学习率为每循环100次，降为上一次的0.01倍，初始为0.01，循环次数epoch为即最低循环256次，batchsize即每次循环内，每次放入计算数据的长度大小，为数据长度的五分之一左右。设定好以上参数后，我们把数据输入模型中，开始建模，并输出预测效果图。
图20
图21
从上图20可以看出，我们建立的优化前LSTM模型预测结果出现了比较明显的“右偏”现象，我们也解释过，如果预测数据出现了相较于原始数据整体往右的现象，说明我们的模型建立是失败的，不能作为预测模型，这很大原因是因为原始数据存在较为明显的自相关性原因。所以我们要对LSTM进行优化，优化的思路跟SVM实质是一样。
4.2：优化后的模型建立
同理SVM模型，我们建立更为复杂的LSTM模型，模型预测结果效果图如下展示：
图21
图22
从上图21观察，同样我们也彻底消除了预测数据结果“右偏”无效的现象，
从上图22的评价指标我们发现，无论是R²值还是MSE值，也是朝着好的方向提升，总体上说明加入因子建模后，模型都得到了很大的提升。
5：两模型的总结
|
模型名称/指标 |
训练集R² |
测试集R² |
总体R² |
训练集MSE |
测试集MSE |
总体MSE |
|
优化前SVM |
0.81 |
0.13 |
0.82 |
3453.45 |
2652.21 |
3290.98 |
|
优化后SVM |
0.86 |
0.43 |
0.87 |
2429.22 |
1753.89 |
2292.28 |
|
优化前LSTM |
0.76 |
0.10 |
0.78 |
4268.02 |
2743.13 |
3985.81 |
|
优化后LSTM |
0.88 |
0.35 |
0.88 |
2210.05 |
1976.67 |
2162.72 |
表1
表1呈现的就是我们模型的主要评价结果，这样很清楚地看出，我们的优化结果提升是非常大的，优化方向和思路是正确客观的。
五：总结
1：总结
本报告主要利用机器学习模型和循环神经网络模型对“订单交易量”时序数据进行预测。首先通过数据清洗和数据预处理过程，得到尽可能“干净无噪声”的数据，再通过描述统计分析和相关检验，得出因子的变化跟交易量的变化的关系强弱，最后再创新地结合相关因子分析，对时序预测模型进行优化，优化过程中涉及的相关指标如R²、MSE值等大大地提升，预测走势更加符合实际。这其中的分析步骤是科学的也是客观的。
2：讨论
本文在处理含有自回归性质的时序数据时，主要是结合相关强因子分析来处理，最终防止预测数据出现“右偏无效”结果，这是一种手段，但也可以通过其它方法，比如对时序数据进行差分处理，最终再还原；再比如建立时间序列相关的线性模型，可以有效防止预测“右偏无效”现象。在选择预测模型的时候，我们也可以建立其它机器学习模型如随机森林RF、Xgboost、多元回归等，深度学习模型如小波神经网络WNN、卷积神经网络CNN等做对比分析，这也是未来我们努力的方向！
六：源码
from sklearn.metrics import mean_squared_error,r2_score
from keras.callbacks import LearningRateScheduler
from keras.models import Sequential
import matplotlib.pyplot as plt
from keras.layers import Dense
from keras.layers import LSTM
from keras import optimizers
import keras.backend as K
import tensorflow as tf
import pandas as pd
import numpy as np
import time
start_time = time.time()
plt.rcParams['font.sans-serif'] = ['SimHei'] ##中文乱码问题！
plt.rcParams['axes.unicode_minus'] = False # 横坐标负号显示问题！
###初始化参数
my_seed = 666 # 随便给出个随机种子
tf.random.set_seed(my_seed) ##运行tf才能真正固定随机种子
path = './'##请正确输入自己文件路径！！！
data_original = pd.read_excel(path + '建模数据.xlsx', index_col=0)
order_data = np.array(data_original['订单量'])
revisedata = np.max(order_data)
order_datanorm = order_data / revisedata
date_ls = data_original['时间']##获取日期列表
num_steps = 5##取预测滑动步长，根据数据量合理选择
test_len = int(len(order_data) * 0.2) ##测试集数量长度
##数据形状转换，很重要！！
def data_format(data, num_steps=3, test_len=5):
print('\033[1;34m{0:*^80}\033[1;0m'.format('数据整合'))
# 根据test_len进行分组
X = np.array([data[i: i + num_steps]
for i in range(len(data) - num_steps)])
y = np.array([data[i + num_steps]
for i in range(len(data) - num_steps)])
train_X, test_X = X[:-test_len], X[-test_len:]
train_y, test_y = y[:-test_len], y[-test_len:]
return train_X, train_y, test_X, test_y
transformer_orderdata = np.reshape(pd.Series(order_datanorm).values, (-1, 1))
train_X, train_y, test_X, test_y = data_format(transformer_orderdata, num_steps, test_len)
print('\033[1;38m原始序列维度信息：%s；转换后训练集X数据维度信息：%s，Y数据维度信息：%s；测试集X数据维度信息：%s，Y数据维度信息：%s\033[0m' % (
transformer_orderdata.shape, train_X.shape, train_y.shape, test_X.shape, test_y.shape))
def buildmylstm(initactivation='relu', ininlr=0.001):
print('\033[1;34m{0:*^80}\033[1;0m'.format('训练模型中。。。'))
nb_lstm_outputs1 = 256 # 隐藏层1神经元个数
nb_lstm_outputs2 = 128 # 隐藏层2神经元个数
nb_time_steps = train_X.shape[1] # 时间序列步长长度
nb_input_vector = train_X.shape[2] # 输入序列特征维度
model = Sequential()
model.add(LSTM(units=nb_lstm_outputs1, input_shape=(nb_time_steps, nb_input_vector), return_sequences=True))
model.add(LSTM(units=nb_lstm_outputs2, input_shape=(nb_time_steps, nb_input_vector)))
model.add(Dense(64, activation=initactivation))
model.add(Dense(32, activation='relu'))
model.add(Dense(test_y.shape[1], activation='tanh'))
lr = ininlr
adam = tf.keras.optimizers.Adam(learning_rate=lr)
def scheduler(epoch): ##编写学习率变化函数
# 每隔epoch，学习率减小为原来的1/10
if epoch % 100 == 0 and epoch != 0:
lr = K.get_value(model.optimizer.lr)
K.set_value(model.optimizer.lr, lr * 0.1)
print('lr changed to {}'.format(lr * 0.1))
return K.get_value(model.optimizer.lr)
model.compile(loss='mse', optimizer=adam, metrics=['mse']) ##根据损失函数性质，回归建模一般选用”距离误差“作为损失函数，分类一般选”交叉熵“损失函数
reduce_lr = LearningRateScheduler(scheduler)
###数据集较少，全参与形式，epochs一般跟batch_size成正比
##callbacks：回调函数，调取reduce_lr
##verbose=0:非冗余打印，即不打印训练过程
batchsize = int(len(order_data) / 5)
epochs = max(128, batchsize * 4) ##最低大循环次数128
model.fit(train_X, train_y, batch_size=batchsize, epochs=epochs, verbose=0, callbacks=[reduce_lr])
return model
def prediction(lstmmodel):
print('\033[1;34m{0:*^80}\033[1;0m'.format('模型评估指标值'))
plt.figure(figsize=(15, 8))
predsinner = lstmmodel.predict(train_X)
predsinner_true = predsinner * revisedata
predsouter = lstmmodel.predict(test_X)
predsouter_true = predsouter * revisedata
init_value = predsinner_true[-1]
plt.plot(order_data, 'r',label='原始值')
Xinner = [i for i in range(num_steps, len(order_data) - test_len)]
plt.plot(Xinner, list(predsinner_true), 'g',label='样本内预测值')
Xouter = [i for i in range(len(order_data) - test_len - 1, len(order_data))]
plt.plot(Xouter, [init_value] + list(predsouter_true), 'b',label='样本外预测值')
set1 = []
set2 = []
L = 10##横坐标刻度个数
for i in np.linspace(1, len(order_data), L):
s = str(date_ls[int(i) - 1])
if len(s) > 10:
set1.append(int(i) - 1)
set2.append(s[:-9])
else:
set1.append(int(i) - 1)
set2.append(s)
plt.xticks(set1, set2, rotation=15, fontsize=8)
plt.legend(fontsize=15)
plt.tick_params(labelsize=15)
plt.xlabel('日期', fontsize=18)
plt.ylabel('订单量', fontsize=18)
plt.title('优化前LSTM模型的预测效果图', fontsize=20)
plt.show()
data_inner = order_data[num_steps:len(predsinner_true) + num_steps]
data_outer = order_data[len(predsinner_true) + num_steps:]
data_all = order_data[num_steps:]
allpredata = list(predsinner_true) + list(predsouter_true)
r2_inner = np.round(r2_score(data_inner, predsinner_true), 2) ##前面是真实数据，后面是预测
r2_outer = np.round(r2_score(data_outer, predsouter_true), 2)
r2_all = np.round(r2_score(data_all, allpredata), 2)
mse_inner = np.round(mean_squared_error(data_inner, predsinner_true), 2) ##MSE函数无所谓顺序
mse_outer = np.round(mean_squared_error(data_outer, predsouter_true), 2)
mse_all = np.round(mean_squared_error(data_all, allpredata), 2)
print('\033[1;35m训练集R方得分：%s，测试集R方得分：%s，总R方得分：%s\033[0m' % (r2_inner, r2_outer, r2_all))
print('\033[1;36m训练集MSE值：%s，测试集MSE值：%s，总MSE值：%s\033[0m' % (mse_inner, mse_outer, mse_all))
return allpredata
mymlstmmodel = buildmylstm()
presult = prediction(mymlstmmodel)
stop_time = time.time()
print('\033[1;34m完整运行结束最终用时（秒）：%s\033[0m' % (round(stop_time - start_time, 3)))
from sklearn.metrics import mean_squared_error,r2_score
from keras.callbacks import LearningRateScheduler
from keras.models import Sequential
import matplotlib.pyplot as plt
from keras.layers import Dense
from keras.layers import LSTM
from keras import initializers
from keras import optimizers
import keras.backend as K
import tensorflow as tf
import pandas as pd
import numpy as np
import time
start_time = time.time()
plt.rcParams['font.sans-serif'] = ['SimHei'] ##中文乱码问题！
plt.rcParams['axes.unicode_minus'] = False # 横坐标负号显示问题！
###初始化参数
my_seed = 666 # 随便给出个随机种子
tf.random.set_seed(my_seed) ##运行tf才能真正固定随机种子
path = './'##请正确输入自己文件路径！！！
data_original = pd.read_excel(path + '建模数据.xlsx', index_col=0)
holiday_work = np.array(data_original['周日_2周六节假日_1班_0'])
order_data = np.array(data_original['订单量'])
meantemp = np.array(data_original['均温度'])
rainy = np.array(data_original['非雨_0雨1'])
holiday_work = holiday_work / np.max(holiday_work)
meantemp = meantemp / np.max(meantemp)
rainy = rainy / np.max(rainy)
revisedata = np.max(order_data)
order_datanorm = order_data / revisedata ##数据规范化
num_steps = 5 ##取序列滑动步长，根据样本实际量来定
test_len = int(len(order_data) * 0.2) ##测试集数量长度
date_ls = data_original['时间']##获取日期列表
##数据形状转换，很重要！！
def data_format(data,ed1,ed2,ed3,ed4,num_steps=3, test_len=10):
print('\033[1;34m{0:*^80}\033[1;0m'.format('数据整合'))
# 根据test_len进行分组
y = np.array([data[i + num_steps]
for i in range(len(data) - num_steps)])
Xed1 = np.array([ed1[i: i + num_steps]
for i in range(len(ed1) - num_steps)])
Xed2 = np.array([ed2[i: i + num_steps]
for i in range(len(ed2) - num_steps)])
Xed3 = np.array([ed3[i: i + num_steps]
for i in range(len(ed3) - num_steps)])
Xed4 = np.array([ed4[i: i + num_steps]
for i in range(len(ed4) - num_steps)])
transformer_inewls_fin = []
for i,j,k,l in zip(Xed1,Xed2,Xed3,Xed4):##对数据整合要小心
inew = list(i) + list(j) + list(k) + list(l)
transformer_inewls = np.reshape(pd.Series(inew).values, (-1, 1))
transformer_inewls_fin.append(transformer_inewls)
transformer_inewls_fin = np.array(transformer_inewls_fin)
train_X, test_X = transformer_inewls_fin[:-test_len], transformer_inewls_fin[-test_len:]
train_y, test_y = y[:-test_len], y[-test_len:]
return train_X, train_y, test_X, test_y
transformer_orderdata = np.reshape(pd.Series(order_datanorm).values, (-1, 1))
train_X, train_y, test_X, test_y = data_format(transformer_orderdata,order_datanorm,meantemp,holiday_work,rainy, num_steps, test_len)
print('\033[1;38m原始序列维度信息：%s，设定步长：%s；转换后训练集X数据维度信息：%s，Y数据维度信息：%s；测试集X数据维度信息：%s，Y数据维度信息：%s\033[0m' % (
transformer_orderdata.shape, num_steps,train_X.shape, train_y.shape, test_X.shape, test_y.shape))
def buildmylstm(initactivation='relu', ininlr=0.001):
print('\033[1;34m{0:*^80}\033[1;0m'.format('训练模型中。。。'))
nb_lstm_outputs1 = 256 # 隐藏层1神经元个数
nb_lstm_outputs2 = 128 # 隐藏层2神经元个数
nb_time_steps = train_X.shape[1] # 时间序列步长长度
nb_input_vector = train_X.shape[2] # 输入序列特征维度
model = Sequential()
model.add(LSTM(units=nb_lstm_outputs1, input_shape=(nb_time_steps, nb_input_vector), return_sequences=True))
model.add(LSTM(units=nb_lstm_outputs2, input_shape=(nb_time_steps, nb_input_vector)))
model.add(Dense(64, activation=initactivation))
model.add(Dense(32, activation='tanh'))
model.add(Dense(test_y.shape[1], activation='tanh'))
lr = ininlr
adam = tf.keras.optimizers.Adam(learning_rate=lr)
def scheduler(epoch): ##编写学习率变化函数
# 每隔epoch，学习率减小为原来的1/10
if epoch % 100 == 0 and epoch != 0:
lr = K.get_value(model.optimizer.lr)
K.set_value(model.optimizer.lr, lr * 0.1)
print('lr changed to {}'.format(lr * 0.1))
return K.get_value(model.optimizer.lr)
model.compile(loss='mse', optimizer=adam, metrics=['mse']) ##根据损失函数性质，回归建模一般选用”距离误差“作为损失函数，分类一般选”交叉熵“损失函数
reduce_lr = LearningRateScheduler(scheduler)
###数据集较少，全参与形式，epochs一般跟batch_size成正比
##callbacks：回调函数，调取reduce_lr
##verbose=0:非冗余打印，即不打印训练过程
batchsize = int(len(order_data) / 5)
epochs = max(128, batchsize * 4) ##最低大循环次数128
model.fit(train_X, train_y, batch_size=batchsize, epochs=epochs, verbose=0, callbacks=[reduce_lr])
return model
def prediction(lstmmodel):
print('\033[1;34m{0:*^80}\033[1;0m'.format('模型评估指标值'))
plt.figure(figsize=(15, 8))
predsinner = lstmmodel.predict(train_X)
predsinner_true = predsinner * revisedata
predsouter = lstmmodel.predict(test_X)
predsouter_true = predsouter * revisedata
init_value = predsinner_true[-1]
plt.plot(order_data, 'r',label='原始值')
Xinner = [i for i in range(num_steps, len(order_data) - test_len)]
plt.plot(Xinner, list(predsinner_true), 'g',label='样本内预测值')
Xouter = [i for i in range(len(order_data) - test_len - 1, len(order_data))]
plt.plot(Xouter, [init_value] + list(predsouter_true), 'b',label='样本外预测值')
set1 = []
set2 = []
L = 10##横坐标刻度个数
for i in np.linspace(1, len(order_data), L):
s = str(date_ls[int(i) - 1])
if len(s) > 10:
set1.append(int(i) - 1)
set2.append(s[:-9])
else:
set1.append(int(i) - 1)
set2.append(s)
plt.xticks(set1, set2, rotation=15, fontsize=8)
plt.legend(fontsize=15)
plt.tick_params(labelsize=15)
plt.xlabel('日期', fontsize=18)
plt.ylabel('订单量', fontsize=18)
plt.title('优化后LSTM模型的预测效果图', fontsize=20)
plt.show()
data_inner = order_data[num_steps:len(predsinner_true) + num_steps]
data_outer = order_data[len(predsinner_true) + num_steps:]
data_all = order_data[num_steps:]
allpredata = list(predsinner_true) + list(predsouter_true)
r2_inner = np.round(r2_score(data_inner, predsinner_true), 2) ##前面是真实数据，后面是预测
r2_outer = np.round(r2_score(data_outer, predsouter_true), 2)
r2_all = np.round(r2_score(data_all, allpredata), 2)
mse_inner = np.round(mean_squared_error(data_inner, predsinner_true), 2) ##MSE函数无所谓顺序
mse_outer = np.round(mean_squared_error(data_outer, predsouter_true), 2)
mse_all = np.round(mean_squared_error(data_all, allpredata), 2)
print('\033[1;35m训练集R方得分：%s，测试集R方得分：%s，总R方得分：%s\033[0m' % (r2_inner, r2_outer, r2_all))
print('\033[1;36m训练集MSE值：%s，测试集MSE值：%s，总MSE值：%s\033[0m' % (mse_inner, mse_outer, mse_all))
return allpredata
mymlstmmodel = buildmylstm()
presult = prediction(mymlstmmodel)
stop_time = time.time()
print('\033[1;34m完整运行结束最终用时（秒）：%s\033[0m' % (round(stop_time - start_time, 3)))
from sklearn.metrics import mean_squared_error,r2_score
import matplotlib.pyplot as plt
from scipy import stats
from sklearn import svm
import pandas as pd
import numpy as np
plt.rcParams['font.sans-serif']=['SimHei']##中文乱码问题！
plt.rcParams['axes.unicode_minus']=False#横坐标负号显示问题！
path = './'##请正确输入自己文件路径！！！
data_original = pd.read_excel(path + '建模数据.xlsx', index_col=0)
holiday_work = np.array(data_original['周日_2周六节假日_1班_0'])
order_data = np.array(data_original['订单量'])
meantemp = np.array(data_original['均温度'])
rainy = np.array(data_original['非雨_0雨1'])
holiday_work = holiday_work / np.max(holiday_work)
meantemp = meantemp / np.max(meantemp)
rainy = rainy / np.max(rainy)
revisedata = 1
order_datanorm = order_data / revisedata ##数据规范化
date_ls = data_original['时间']##获取日期列表
def svmreg_model(X1,Y1):
svrmodel = svm.SVR(kernel='linear',C=10)##多项式核函数
# c_can = np.linspace(0.001,0.01,10)
# degree_can = np.linspace(1,4,10)
# svrmodel = GridSearchCV(svrmodel,param_grid={'C':c_can,'degree':degree_can},cv=10)
svrmodel.fit(X1,Y1)
return svrmodel
##数据形状转换
def data_format(data,tl,ns=10):
print('\033[1;34m{0:*^80}\033[1;0m'.format('数据划分'))
# 根据test_len进行分组
X = np.array([data[i: i + ns] for i in range(len(data) - ns)])
y = np.array([data[i + ns] for i in range(len(data) - ns)])
train_X, test_X = X[:-tl], X[-tl:]
train_y, test_y = y[:-tl], y[-tl:]
print('\033[1;33m数据总长度：%s，当设定步数为%s步时，训练集长度：%s，测试集长度：%s\033[0m'%(len(data),ns,len(train_y),len(test_y)))
return train_X,train_y,test_X,test_y,X
num_steps = 5 ##取序列滑动步长，根据样本实际量来定
test_len = int(len(order_data) * 0.25) ##测试集数量长度
rescopy = order_data.copy()
# S = pd.Series(order_data).std()
# m = pd.Series(order_data).mean()
# res2 = (np.array(res2) - m) / S#数据规范化
res3 = data_format(order_data,test_len,num_steps)
res4 = svmreg_model(res3[0],res3[1])
def evaluatemodel(data):
print('\033[1;34m{0:*^80}\033[1;0m'.format('模型评估指标值'))
preinner = res4.predict(res3[0]) * revisedata
preouter = res4.predict(res3[2]) * revisedata
init_value = preinner[-1]
plt.figure(figsize=(15, 8))
Xinner = [i for i in range(num_steps, len(data) - test_len)]
Xouter = [i for i in range(len(data) - test_len - 1, len(data))]
allpredata = list(preinner) + list(preouter)
plt.plot(data, 'r',label='预处理后原始值')
plt.plot(Xinner, preinner, 'g',label='样本内预测值')
plt.plot(Xouter, [init_value] + list(preouter), 'b',label='样本外预测值')
plt.title('优化前SVM模型的预测效果图',fontsize=20)
plt.legend(fontsize=18)
plt.tick_params(labelsize=15)
set1 = []
set2 = []
if len(data) <= 16:
L = len(data)
else:
L = 10##横坐标刻度个数
for i in np.linspace(1, len(data), L):
s = str(date_ls[int(i) - 1])
if len(s) > 10:
set1.append(int(i) - 1)
set2.append(s[:-9])
else:
set1.append(int(i) - 1)
set2.append(s)
plt.xticks(set1,set2,rotation=15,fontsize=8)
plt.legend(fontsize=15)
plt.tick_params(labelsize=15)
plt.xlabel('日期',fontsize=18)
plt.ylabel('订单量',fontsize=18)
plt.show()
data_inner = data[num_steps:len(preinner) + num_steps]
data_outer = data[len(preinner) + num_steps:]
data_all = data[num_steps:]
r2_inner = np.round(r2_score(data_inner,preinner),2)##前面是真实数据，后面是预测
r2_outer = np.round(r2_score(data_outer,preouter),2)
r2_all = np.round(r2_score(data_all,allpredata),2)
mse_inner = np.round(mean_squared_error(data_inner,preinner),2)##MSE函数无所谓顺序
mse_outer = np.round(mean_squared_error(data_outer,preouter),2)
mse_all = np.round(mean_squared_error(data_all,allpredata),2)
print('\033[1;35m训练集R方得分：%s，测试集R方得分：%s，总R方得分：%s\033[0m'%(r2_inner,r2_outer,r2_all))
print('\033[1;36m训练集MSE值：%s，测试集MSE值：%s，总MSE值：%s\033[0m'%(mse_inner,mse_outer,mse_all))
evaluatemodel(rescopy)
from sklearn.metrics import mean_squared_error,r2_score
import matplotlib.pyplot as plt
from scipy import stats
from sklearn import svm
import pandas as pd
import numpy as np
plt.rcParams['font.sans-serif']=['SimHei']##中文乱码问题！
plt.rcParams['axes.unicode_minus']=False#横坐标负号显示问题！
path = './'##请正确输入自己文件路径！！！
data_original = pd.read_excel(path + '建模数据.xlsx', index_col=0)
holiday_work = np.array(data_original['周日_2周六节假日_1班_0'])
order_data = np.array(data_original['订单量'])
meantemp = np.array(data_original['均温度'])
rainy = np.array(data_original['非雨_0雨1'])
holiday_work = holiday_work / 1
meantemp = meantemp / 1
rainy = rainy / 1
revisedata = 1
order_datanorm = order_data / revisedata ##数据规范化
date_ls = data_original.index.to_list()##获取日期列表
def svmreg_model(X1,Y1):
svrmodel = svm.SVR(kernel='linear',C=10)##多项式核函数
# c_can = np.linspace(0.001,0.01,10)
# degree_can = np.linspace(1,4,10)
# svrmodel = GridSearchCV(svrmodel,param_grid={'C':c_can,'degree':degree_can},cv=10)
svrmodel.fit(X1,Y1)
return svrmodel
##数据形状转换
def data_format(data,ed1,ed2,ed3,ed4,num_steps=3,test_len=10):
print('\033[1;34m{0:*^80}\033[1;0m'.format('数据划分'))
y = np.array([data[i + num_steps]
for i in range(len(data) - num_steps)])
Xed1 = np.array([ed1[i: i + num_steps]
for i in range(len(ed1) - num_steps)])
Xed2 = np.array([ed2[i: i + num_steps]
for i in range(len(ed2) - num_steps)])
Xed3 = np.array([ed3[i: i + num_steps]
for i in range(len(ed3) - num_steps)])
Xed4 = np.array([ed4[i: i + num_steps]
for i in range(len(ed4) - num_steps)])
transformer_inewls_fin = []
for i,j,k,l in zip(Xed1,Xed2,Xed3,Xed4):##对数据整合要小心
inew = list(i) + list(j) + list(k) + list(l)
transformer_inewls_fin.append(inew)
transformer_inewls_fin = np.array(transformer_inewls_fin)
train_X, test_X = transformer_inewls_fin[:-test_len], transformer_inewls_fin[-test_len:]
train_y, test_y = y[:-test_len], y[-test_len:]
return train_X, train_y, test_X, test_y
num_steps = 5 ##取序列滑动步长，根据样本实际量来定
test_len = int(len(order_data) * 0.25) ##测试集数量长度
rescopy = order_data.copy()
# S = pd.Series(order_data).std()
# m = pd.Series(order_data).mean()
# res2 = (np.array(res2) - m) / S#数据规范化
train_X, train_y, test_X, test_y = data_format(order_data,order_datanorm,meantemp,holiday_work,rainy, num_steps, test_len)
print('\033[1;38m原始序列维度信息：%s，设定步长：%s；转换后训练集X数据维度信息：%s，Y数据维度信息：%s；测试集X数据维度信息：%s，Y数据维度信息：%s\033[0m' % (
data_original.shape, num_steps,train_X.shape, train_y.shape, test_X.shape, test_y.shape))
res4 = svmreg_model(train_X,train_y)
def evaluatemodel(data):
print('\033[1;34m{0:*^80}\033[1;0m'.format('模型评估指标值'))
preinner = res4.predict(train_X) * revisedata
preouter = res4.predict(test_X) * revisedata
init_value = preinner[-1]
plt.figure(figsize=(15, 8))
Xinner = [i for i in range(num_steps, len(data) - test_len)]
Xouter = [i for i in range(len(data) - test_len - 1, len(data))]
allpredata = list(preinner) + list(preouter)
plt.plot(data, 'r',label='预处理后原始值')
plt.plot(Xinner, preinner, 'g',label='样本内预测值')
plt.plot(Xouter, [init_value] + list(preouter), 'b',label='样本外预测值')
plt.legend(fontsize=18)
plt.tick_params(labelsize=15)
set1 = []
set2 = []
if len(data) <= 16:
L = len(data)
else:
L = 10##横坐标刻度个数
for i in np.linspace(1, len(data), L):
s = str(date_ls[int(i) - 1])
if len(s) > 10:
set1.append(int(i) - 1)
set2.append(s[:-9])
else:
set1.append(int(i) - 1)
set2.append(s)
plt.xticks(set1,set2,rotation=15,fontsize=8)
plt.legend(fontsize=15)
plt.title('优化后SVM模型的预测效果图', fontsize=20)
plt.tick_params(labelsize=15)
plt.xlabel('日期',fontsize=18)
plt.ylabel('订单量',fontsize=18)
plt.show()
data_inner = data[num_steps:len(preinner) + num_steps]
data_outer = data[len(preinner) + num_steps:]
data_all = data[num_steps:]
r2_inner = np.round(r2_score(data_inner,preinner),2)##前面是真实数据，后面是预测
r2_outer = np.round(r2_score(data_outer,preouter),2)
r2_all = np.round(r2_score(data_all,allpredata),2)
mse_inner = np.round(mean_squared_error(data_inner,preinner),2)##MSE函数无所谓顺序
mse_outer = np.round(mean_squared_error(data_outer,preouter),2)
mse_all = np.round(mean_squared_error(data_all,allpredata),2)
print('\033[1;35m训练集R方得分：%s，测试集R方得分：%s，总R方得分：%s\033[0m'%(r2_inner,r2_outer,r2_all))
print('\033[1;36m训练集MSE值：%s，测试集MSE值：%s，总MSE值：%s\033[0m'%(mse_inner,mse_outer,mse_all))
evaluatemodel(rescopy)
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns
import pandas as pd
import numpy as np
import time
start_time = time.time()
plt.rcParams['font.sans-serif']=['SimHei']##中文乱码问题！
plt.rcParams['axes.unicode_minus']=False#横坐标负号显示问题！
path = './'##请正确输入自己文件路径！！！
data = pd.read_excel(path + '某地区某年共享电动车需求量数据.xlsx')
print('\033[1;34m{0:*^80}\033[1;0m'.format('缺失值查询'))
datashape = data.shape
print('\033[1;33m原数据形状：\033[0m',datashape)
col_ls = list(data)[1:]
originallen = datashape[0]
index_ls = data['时间']
def lookup_missvalue(colls):
print('\033[1;34m{0:*^80}\033[1;0m'.format('缺失值查询'))
missvalue_ls = []
datashape = data.shape
print('\033[1;33m原数据形状：\033[0m',datashape)
for i in colls:##性别和编号不做计算
datatemp = data[i].dropna()
if len(datatemp) < datashape[0]:
print('\033[1;31m特征“%s”的数据长度为：%s，小于%s，为缺失数据特征\033[0m' % (i, len(datatemp),datashape[0]))
missvalue_ls.append(i)
else:
print('\033[1;32m特征“%s”的数据长度为：%s，等于%s，正常特征\033[0m' % (i, len(datatemp), datashape[0]))
return missvalue_ls
res = lookup_missvalue(col_ls)
def fill_missvalue(misscol_ls):
print('\033[1;34m{0:*^80}\033[1;0m'.format('缺失值补充'))
for i in misscol_ls:
med = data[i].median()
data[i] = data[i].fillna(med) ##缺失值用中位数填充
print('\033[1;35m缺失值特征“%s”数据用中位数%s替换\033[0m' % (i,med))
print(data)
fill_missvalue(res)
def substitute_outliner(data,original_ls,yl,tit,mul1,mul2,num_ls):
print('\033[1;34m{0:*^80}\033[1;0m'.format('特征数据“%s”异常值处理'%yl))
copydata = original_ls.copy()
for i in data[5]:##获取异常值位置
del original_ls[i]##按照位置删除最最保险
upvalue = data[2][i]
downvalue = data[3][i]
originalvalue = copydata[i]
if originalvalue <= downvalue:
changevaluetemp = downvalue * 0.5 + upvalue * 0.25
if changevaluetemp >= upvalue:
changevalue = upvalue
else:
changevalue = max(changevaluetemp,downvalue)##最小不低于downvalue，防止upvalue过小
original_ls.insert(i,round(changevalue,1))
print('\033[1;31m异常数据编号为（原始值%s）：%s，替换的修正值为：%s，downvalue：%s\033[0m' % (originalvalue, str(num_ls[i])[:-9], round(changevalue,2),round(downvalue,2)))
if originalvalue >= upvalue:
changevaluetemp = upvalue * 0.5 + downvalue * 0.25
if changevaluetemp <= downvalue:
changevalue = downvalue
else:
changevalue = min(changevaluetemp,upvalue)##最大不超过upvalue，防止downvalue过大
original_ls.insert(i,round(changevalue,1))
print('\033[1;31m异常数据编号为（原始值%s）：%s，替换的修正值为：%s，upvalue：%s\033[0m'%(originalvalue,str(num_ls[i])[:-9],round(changevalue,2),round(upvalue,2)))
plt.figure(figsize=(15, 8))
plt.plot(list(copydata), 'purple', label='原始数值')
plt.plot(data[1], 'grey', label='拟合曲线')
plt.plot(data[2], 'r', label='拟合上限(%ssigm)曲线' % mul1)
plt.plot(data[3], 'b', label='拟合下限(%ssigm)曲线' % mul2)
plt.plot(original_ls, 'g', label='处理完后的数据')
tuple1 = [] ##定义两个刻度值"元组"
tuple2 = []
L = 10 ##给出刻度数据量
for i in np.linspace(1, len(copydata), L):
s = str(num_ls[int(i) - 1])
if len(s) > 10:
tuple1.append(int(i) - 1)
tuple2.append(s[:-9])
else:
tuple1.append(int(i) - 1)
tuple2.append(s)
plt.xticks(tuple1, tuple2, rotation=15, fontsize=8)
plt.legend(fontsize=15)
plt.tick_params(labelsize=15)
plt.xlabel('日期', fontsize=18)
plt.ylabel(yl, fontsize=18)
plt.title(tit, fontsize=20)
plt.show()
return original_ls
def check_outlier(value,mul1,mul2,yl):
print('\033[1;34m{0:*^80}\033[1;0m'.format('特征数据“%s”预处理'%yl))
x = [j for j in range(len(value))]
coeffs = np.polyfit(x, value, 5) ##专门求多项式估计参数的函数，根据实际数据波动取相应的阶数
p = np.poly1d(coeffs) # 一元估计参数
sigm = p(x).std()
sigm_up = p(x) + mul1 * sigm
sigm_down = p(x) - mul2 * sigm
outliner = value[(value < sigm_down) | (value > sigm_up)]#条件筛选
print('\033[1;31m原数据长度：%s，异常数据：%s\033[0m'%(len(value),len(outliner)))
print(outliner)
return list(outliner),p(x),sigm_up,sigm_down,x,list(outliner.index)
for i in ['订单量','均温度']:##性别和编号不做计算
mul1 = 3
mul2 = 3##由于数据很大是人为的，这里mul1和mul2参数慎改
res1 = check_outlier(data[i],mul1,mul2,i)
res2 = substitute_outliner(res1,list(data[i]),i,'“%s”数据预处理图'%i,mul1,mul2,index_ls)
data[i] = res2
data.to_excel(path + '建模数据.xlsx')
##对于周末、周六（节假日）、工作日进行作图可视化
numsunday = data[data.周日_2周六节假日_1班_0==2]['订单量']
numsat_holiday = data[data.周日_2周六节假日_1班_0==1]['订单量']
numwork = data[data.周日_2周六节假日_1班_0==0]['订单量']
numsundayt = round(numsunday.mean(),0)
numsat_holidayt = round(numsat_holiday.mean(),0)
numworkt = round(numwork.mean(),0)
plt.figure(figsize=(15, 8))
num_ls1 = [numsundayt,numsat_holidayt,numworkt]
label_ls1 = ['周日','周六和法定节假日','工作日']
for i,j,l in zip(num_ls1,range(1,len(num_ls1) + 1),label_ls1):
plt.bar(j,i,width=0.5,label=l)
plt.xticks((1,2,3),('周日','周六和法定节假日','工作日'),rotation=15,fontsize=18)
plt.ylabel('平均每天需求量',fontsize=18)
plt.tick_params(labelsize=18)
plt.legend(fontsize=15)
plt.grid(axis='y')
##对于雨天、非雨天进行作图可视化
plt.figure(figsize=(15, 8))
numrainy = data[data.非雨_0雨1==1]['订单量']
numnorainy = data[data.非雨_0雨1==0]['订单量']
numrainyt = round(numrainy.mean(),0)
numnorainyt = round(numnorainy.mean(),0)
num_ls2 = [numrainyt,numnorainyt]
label_ls2 = ['雨天','非雨天']
for i,j,l in zip(num_ls2,range(1,len(num_ls2) + 1),label_ls2):
plt.bar(j,i,width=0.75,label=l)
plt.xticks((1,2),('雨天','非雨天'),rotation=15,fontsize=18)
plt.ylabel('平均每天需求量',fontsize=18)
plt.tick_params(labelsize=18)
plt.legend(fontsize=15)
plt.grid(axis='y')
##对于气温进行作图可视化
plt.figure(figsize=(15, 8))
meant1 = data[data.均温度<-0]['订单量']
meant2 = data[(data.均温度>=0) & (data.均温度<10)]['订单量']
meant3 = data[(data.均温度>=10) & (data.均温度<20)]['订单量']
meant4 = data[(data.均温度>=20) & (data.均温度<30)]['订单量']
meant5 = data[data.均温度>30]['订单量']
numtemp1 = round(meant1.mean(),0)
numtemp2 = round(meant2.mean(),0)
numtemp3 = round(meant3.mean(),0)
numtemp4 = round(meant4,0)
numtemp5 = round(meant5,0)
num_ls3 = [numtemp1,numtemp2,numtemp3,numtemp4,numtemp5]
label_ls3 = ['<0度','[0,10)度','[10,20)度','[20,30)度','>30度']
for i,j,l in zip(num_ls3,range(1,len(num_ls3) + 1),label_ls3):
plt.bar(j,i,width=0.75,label=l)
plt.ylabel('平均每天需求量',fontsize=18)
plt.tick_params(labelsize=18)
plt.legend(fontsize=15)
plt.grid(axis='y')
plt.show()
plt.figure(figsize=(15, 8)) ##年龄分布图统计
sns.distplot(data['订单量'], rug=False, hist=True)
plt.tick_params(labelsize=18)
plt.ylabel('密度', fontsize=18)
plt.xlabel('订单交易量', fontsize=18)
# plt.legend(fontsize=15)
plt.show()
print(data['订单量'].describe())
####数据检验
test_meant = stats.kruskal(meant1,meant2,meant3,meant4,meant5)
print('\033[1;32m“温度因素”是否对订单量产生影响的检验p值：%s\033[0m'%test_meant[1])
test_day = stats.kruskal(numsunday,numsat_holiday,numwork)
print('\033[1;32m“是否假期因素”是否对订单量产生影响的检验p值：%s\033[0m'%test_day[1])
test_rain = stats.kruskal(numrainy,numnorainy)
print('\033[1;32m“是否雨天因素”是否对订单量产生影响的检验p值：%s\033[0m'%test_rain[1])
import statsmodels.api as sm
def testwhitenoise(data):
m = 15# 检验10个自相关系数
acf,q,p = sm.tsa.acf(data,nlags=m,qstat=True)
out = np.c_[range(1,m+1),acf[1:],q,p]
output = pd.DataFrame(out,columns=['lag','自相关系数','统计量Q值','p_values'])
output = output.set_index('lag')# 设置第一列索引名称,可省略重复索引列1
print(output)
testwhitenoise(data['订单量'])
from statsmodels.tsa.seasonal import seasonal_decompose
ts_decomposition = seasonal_decompose(data['订单量'], period=14)##数据分解，一般period设为12，一个完整年月份
plt.figure(figsize=(15, 8))
ts_decomposition.plot()
plt.xlabel('数据结构分解图', fontsize=10)
plt.show()
标签：电动车,plt,报告,len,test,np,033,共享,data From： https://www.cnblogs.com/ctcqq/p/16997167.html