python实训——回归类型数据挖掘任务

1. 回归类型数据挖掘任务

基于ARIMA和多层神经网络模型的地铁站点日客流量预测。有郑州市2015年8月-11月各地铁闸机刷卡数据集。对每日各地铁站的客流量进行分析并进行可视化。基于上一步的分析结果，分别采用ARIMA模型和多层神经网络模型对数据进行建模，训练优化模型并分别给出评估指标。

原始数据如下所示，共分为4个文件，近15GB。

2.1总体流程

2.2数据分析

2.2.1读取数据

使用pandas读取下列文件，并将读取到的DataFrame进行合并。

acc_08_final_mini.csv

acc_09_final_mini.csv

acc_10_final_mini.csv

acc_11_final_mini.csv

得到的DataFrame如图所示：

代码：

import pandas as pd
# 创建一个空的DataFrame列表来存储所有读取的DataFrame  
df_list = []  
  
# 循环读取所有CSV文件  
for year in range(8, 12):  # 从2008年到2011年  
    filename = fr'C:\Users\zjl15\PycharmProjects\pythonProject1\shixundata\2\acc_{year:02d}_final_mini.csv'  # 注意fr的使用  
    df = pd.read_csv(filename)  # 读取CSV文件  
    df_list.append(df)  # 将读取的DataFrame添加到列表中  
  
# 使用concat函数进行纵向合并（堆叠）  
combined_df = pd.concat(df_list, ignore_index=True)  
  
# 显示合并后的DataFrame  
print(combined_df)  
  

特征含义：

特征名	TRADE_TYPE	TRADE_ADDRESS	TRADE_DATE
含义	交易类型 21：进站 22：出站	交易站点	交易时间

2.2.2分析数据

2.2.2.1日期处理

提取TRADE_DATE中的日期部分作为新的TRADE_DATE的内容。

建议使用pandas库的apply函数。

2.2.2.2分组统计

统计各个站点每天的客流量，并将客流量数据作为1列加入到原数据中，列名为COUNT。建议使用pandas库的groupby函数和transform函数。

2.2.2.3删除TRADE_TYPE列

因本次分析不考虑交易类型这个因素，因此要删掉TRADE_TYPE列。建议使用pandas库的drop函数。

2.2.2.4排序去重

将数据按照TRADE_ADDRESS和TRADE_DATE两列进行排序，并进行去重处理。建议使用pandas库的sort_values函数和drop_duplicates函数。

# 将 TRADE_DATE 列转换为 datetime 类型
combined_df['TRADE_DATE'] = pd.to_datetime(combined_df['TRADE_DATE'], format='%Y-%m-%d-%H.%M.%S.%f')
# 使用 apply 函数提取日期部分并创建新的列
combined_df['DATE_ONLY'] = combined_df['TRADE_DATE'].apply(lambda x: x.date())
#去除TRADE_DATE列
df=combined_df.drop(columns=['TRADE_DATE'])
# 使用 groupby 和 transform 计算每个站点每天的交易数量并将结果作为新列加入原数据
df['COUNT'] = df.groupby(['DATE_ONLY', 'TRADE_ADDRESS'])['TRADE_TYPE'].transform('count')
#去除TRADE_TYPE列
df=df.drop(columns=['TRADE_TYPE'])
# 按照 TRADE_ADDRESS 和 DATE_ONLY 进行排序
df_sorted = df.sort_values(by=['TRADE_ADDRESS', 'DATE_ONLY'])
# 去重处理
df_unique = df_sorted.drop_duplicates(subset=['TRADE_ADDRESS', 'DATE_ONLY'])
# 根据日期和站点对数据进行分组，并计算客流量总和
grouped_data = df_unique.groupby(['DATE_ONLY', 'TRADE_ADDRESS'])['COUNT'].sum().unstack()

 

2.2.2.5绘制所有站点8月1日-11月30日的客流量折线图

# 绘制折线图
grouped_data.plot(figsize=(12, 8))
plt.title('Daily Traffic for Each Station')
plt.xlabel('Date')
plt.ylabel('Traffic Count')
plt.legend(title='Station', bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

 

2.3模型构建

2.3.1 ARIMA模型（以121站点为例）

以2.2.2.4步骤中的数据为基础，筛选出TRADE_ADDRESS=121的行。

以121站点客流量数据为例，使用时间序列分析方法进行建模。

# 从原始数据中筛选出TRADE_ADDRESS=121的行
data_121 = df_unique[df_unique['TRADE_ADDRESS'] == 121]
#重新定义一个变量，将数据备份
data = data_121
#删除多余的列
data_121 = data_121.drop(columns=['TRADE_ADDRESS','DATE_ONLY'])

2.3.1.1划分训练集和测试集

共122条记录，前115条为训练集，后7条为测试集

# 划分训练集和测试集
train_set = data_121.iloc[:115]
test_set = data_121.iloc[115:]

2.3.1.2平稳性检验

对训练集进行平稳性检验，分别采用时序图、自相关图、ADF单位根检测方法。

（1）绘制时序图

（2）绘制自相关图

（3）ADF单位根检测

通过检验，发现该序列为非平稳序列。如想使用ARIMA模型进行建模分析，需处理转换为平稳序列。

# 绘制时序图
plt.figure(figsize=(12, 6))
plt.plot(train_set['COUNT'])
plt.title('Traffic Count at Station 121 (Training Set)')
plt.xlabel('Date')
plt.ylabel('Traffic Count')
plt.show()
# 绘制自相关图
fig, ax = plt.subplots(2, 1, figsize=(12, 8))
sm.graphics.tsa.plot_acf(train_set['COUNT'], lags=30, ax=ax[0])
sm.graphics.tsa.plot_pacf(train_set['COUNT'], lags=30, ax=ax[1])
plt.show()
# ADF单位根检测
result = adfuller(train_set['COUNT'])
print('ADF Statistic:', result[0])
print('p-value:', result[1])
print('Critical Values:', result[4])

 

 

2.3.1.3差分处理

使用2阶差分对训练集进行处理。

对差分处理后的训练集再次进行平稳性检验，分别采用时序图、自相关图、ADF单位根检测方法。

# 对训练集进行2阶差分处理
train_set_diff = train_set['COUNT'].diff().diff().dropna()



# 绘制2阶差分后的时序图
plt.figure(figsize=(12, 6))
plt.plot(train_set_diff)
plt.title('Differenced Traffic Count at Station 121 (Training Set)')
plt.xlabel('Date')
plt.ylabel('Differenced Traffic Count')
plt.show()
# 绘制2阶差分后的自相关图
fig, ax = plt.subplots(2, 1, figsize=(12, 8))
sm.graphics.tsa.plot_acf(train_set_diff, lags=30, ax=ax[0])
sm.graphics.tsa.plot_pacf(train_set_diff, lags=30, ax=ax[1])
plt.show()
# ADF单位根检测（2阶差分后）
result_diff = adfuller(train_set_diff)
print('ADF Statistic (after differencing):', result_diff[0])
print('p-value (after differencing):', result_diff[1])
print('Critical Values (after differencing):', result_diff[4])

 

 

2.3.1.4纯随机性检验

对差分处理后的训练集进行纯随机性检验。

# 计算自相关系数
acf = sm.tsa.acf(train_set_diff, fft=False)
# 绘制自相关函数图
plt.figure(figsize=(12, 6))
plt.bar(range(len(acf)), acf)
plt.xlabel('Lag')
plt.ylabel('ACF')
plt.title('Autocorrelation Function')
plt.show()

 

2.3.1.5 ARIMA模型建模

选择合适的参数q和q来构建ARIMA模型。

2.3.1.6 ARIMA模型训练

1. 使用差分处理后的数据对ARIMA模型进行训练
2. 绘制训练集真实值和预测值对比图

参考代码：

# 定义 ARIMA 模型的参数 (p, d, q)
p = 6  # 自回归阶数
d = 2  # 差分阶数
q = 6  # 移动平均阶数
# 训练 ARIMA 模型
model = ARIMA(train_set, order=(p, d, q))
model_fit = model.fit()
# 使用训练好的模型进行预测

# 绘制训练集真实值和预测值对比图
plt.figure(figsize=(12, 6))
plt.plot(train_set, label='Actual')
plt.plot(forecast, label='Predicted', color='red')
plt.title('Comparison of Actual and Predicted Values')
plt.xlabel('Time')
plt.ylabel('Value')
plt.legend()
plt.show()

2.3.2多层神经网络模型

2.3.2.1 特征构造

构造多层神经网络模型所需特征。

1. 筛选121站点数据

以2.2.2.4步骤中的数据为基础，筛选出TRADE_ADDRESS=121的行。

1. 删掉TRADE_ADDRESS列

因所有行均为121站点的数据，所以TRADE_ADDRESS没有意义，删掉。建议使用pandas库的drop函数。

1. 从TRADE_DATE列中提取出星期几数据

对TRADE_DATE列进行计算，根据日期推算出星期几，然后赋值给新的一列TRADE_DAY。建议使用pandas库的apply函数。

1. 从TRADE_DAY列中提取出周末数据

对TRADE_DAY列进行计算，得出是否为周末，然后赋值给新的一列WEEKEND（1：周末，0：工作日）。建议使用pandas库的apply函数。

1. 从TRADE_DATE列中提取出月份数据

对TRADE_DATE列进行计算，提取出月份信息，然后赋值给新的一列MONTH。建议使用pandas库的apply函数。

1. 从TRADE_DATE列中提取出日数据

对TRADE_DATE列进行计算，提取出日信息，然后赋值给新的一列DAY。建议使用pandas库的apply函数。

1. 删掉TRADE_DATE列

TRADE_DATE里的信息已经提取完毕，可以删掉。建议使用pandas库的drop函数。

1. 对数据的列进行重新排序

顺序为：['MONTH','DAY','TRADE_DAY','WEEKEND', 'COUNT']

#删掉TRADE_ADDRESS列
df = data.drop(columns=['TRADE_ADDRESS'])
# 将 DATE_ONLY 列转换为日期时间类型
df['DATE_ONLY'] = pd.to_datetime(df['DATE_ONLY'])
# 使用 apply 函数计算星期几，并赋值给新列 DAY_OF_WEEK
df['DAY_OF_WEEK'] = df['DATE_ONLY'].dt.day_name()
# 使用 dt.dayofweek + 1 将星期几转换为数字表示，并赋值给新列 DAY_OF_WEEK_NUM
df['DAY_OF_WEEK_NUM'] = df['DATE_ONLY'].dt.dayofweek + 1
#删除DAY_OF_WEEK
df2 = df.drop(columns=['DAY_OF_WEEK'])
# 将 TRADE_DATE 列转换为日期时间类型
df2['DATE_ONLY'] = pd.to_datetime(df2['DATE_ONLY'])
# 使用 apply 函数结合 dt.month 提取出月份信息，并赋值给新列 MONTH
df2['MONTH'] = df2['DATE_ONLY'].apply(lambda x: x.month)
# 使用 apply 函数结合 dt.day 提取出日信息，并赋值给新列 DAY
df2['DAY'] = df2['DATE_ONLY'].apply(lambda x: x.day)
#删除DATE_ONLY
df2 = df2.drop(columns=['DATE_ONLY'])
# 自定义函数来判断是否为周末
def is_weekend(day_of_week):
    return 1 if day_of_week >= 6 else 0
# 使用 apply 函数调用自定义函数对 DAY_OF_WEEK_NUM 列进行计算，生成 WEEKEND 列
df2['WEEKEND'] = df2['DAY_OF_WEEK_NUM'].apply(is_weekend)
#按指定顺序排好['MONTH', 'DAY', 'DAY_OF_WEEK_NUM', 'WEEKEND', 'COUNT']
df2 = df2[['MONTH', 'DAY', 'DAY_OF_WEEK_NUM', 'WEEKEND', 'COUNT']]

1. 绘制数据的多变量联合分布图

使用seaborn库的pairplot方法绘制数据的联合分布图。

# 绘制多变量联合分布图
sns.pairplot(df2)
plt.show()

 

2.3.2.2划分训练集和测试集

1. 提取特征x和目标值y

根据2.3.2.1构建出的数据集，提取出特征x和目标值y。

1. 划分训练集和测试集

共122条记录，前115条为训练集，后7条为测试集。

# 提取特征 x
X = df2[['MONTH', 'DAY', 'DAY_OF_WEEK_NUM', 'WEEKEND']]  # 选择需要作为特征的列
# 提取目标值 y
y = df2['COUNT']  # 选择目标列
# 划分训练集和测试集
# 划分训练集和测试集
X_train = X[:115]  # 前115条作为训练集
X_test = X[115:]   # 后7条作为测试集
y_train = y[:115]  # 对应的目标值作为训练集
y_test = y[115:]   # 对应的目标值作为测试集

2.3.2.3模型构建

模型参数：

参考代码

#构建多层神经网络模型
model = tf.keras.Sequential([
    tf.keras.layers.Dense(64,activation='relu',input_shape=[len(X_train.keys())]),
    tf.keras.layers.Dropout(0.3),
    tf.keras.layers.Dense(64,activation='relu'),
    tf.keras.layers.Dropout(0.3),
    tf.keras.layers.Dense(64,activation='relu'),
    tf.keras.layers.Dropout(0.3),
    tf.keras.layers.Dense(1)
])
# 编译模型
model.compile(optimizer='adam', loss='mse', metrics=['mae','mse'])
model.summary()  
# 训练模型并保存历史数据
history = model.fit(X_train, y_train, epochs=3000, validation_split=0.3)

 

2.3.2.4模型编译

编译参数参考：

# 优化器

optimizer='adam'

# 损失函数

loss='mse'

# 评估指标

metrics=['mae','mse']

2.3.2.5模型训练

使用fit函数对训练集进行拟合训练，并将训练过程中产生的历史数据history保存至变量中。

训练参数参考：

# 迭代次数

epochs=1000

# 验证集比例

validation_split=0.2

2.3.2.6训练过程可视化

对history中保存下来的训练过程中的mae和mse的变化情况进行绘图。

参考代码：

# 获取训练过程中的指标数值
mae = history.history['mae']
mse = history.history['loss']
val_mae = history.history['val_mae']
val_mse = history.history['val_mse']
#可视化mae、mse
plt.plot(mae, label='MAE')
plt.plot(val_mae, label='Validation MAE')
plt.title('MAE vs')
plt.legend()
plt.show()
plt.plot(mse, label='MSE')
plt.plot(val_mse, label='Validation MSE')
plt.title('MSE')
plt.legend()
plt.show()
print('MAE:', mae)
print('MSE:', mse)

 

2.4模型评估

2.4.1 ARIMA模型评估

2.4.1.1测试

1. 绘制测试集真实值和预测值对比图
2. 计算模型在测试集上的MAE,MSE,RMSE,MAPE

# 获取训练过程中的指标数值
mae = history.history['mae']
mse = history.history['loss']
val_mae = history.history['val_mae']
val_mse = history.history['val_mse']

plt.plot(mae, label='MAE')
plt.plot(val_mae, label='Validation MAE')
plt.title('MAE vs')
plt.legend()
plt.show()
plt.plot(mse, label='MSE')
plt.plot(val_mse, label='Validation MSE')
plt.title('MSE')
plt.legend()
plt.show()

# history.plot()
# plt.plot(history)
# plt.show()
print('MAE:', mae)
print('MSE:', mse)


y_arima = model_fit.forecast(7)
print('预测未来'+str(len(y_arima))+'天的客流量:\n',y_arima)
y_arima.index = test_set.index
# # # 创建一个新的DataFrame来保存预测值，索引与测试集相同  
# test_set = pd.read_csv(r'test.csv')
# print(test_set)
# test_set.to_csv('test_set.xlsx', index=False)
# y_arima = pd.Series(y_arima, index=test_set.index)  
test_set['p'] = y_arima
plt.figure()  
plt.plot(test_set, label='Actual')  
# plt.plot(y_arima, label='Predicted', color='red')  

plt.legend(['y_true','y_arima'])  
plt.show()

 

 

2.4.2多层神经网络模型评估

2.4.2.1测试

1. 根据测试集的预测结果和真实值、ARIMA模型预测结果进行对比并绘图

参考代码：

y_test = pd.DataFrame(y_test)
y_mlp = pd.DataFrame(y_mlp)
y_mlp.index = y_arima.index
y_test['P1']=y_arima
y_test['P2']=y_mlp
print(y_test)
plt.figure()  
plt.plot(y_test, label='Actual')  
plt.legend(['y_true','y_arima','y_mlp'])  
plt.show()

 

 

完整代码：

import pandas as pd  
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.vector_ar.var_model import VAR
from sklearn import metrics
import numpy as np
import tensorflow as tf
from statsmodels.tsa.arima.model import ARIMA
import matplotlib.pyplot as plt
import seaborn as sns
#%%
import pandas as pd
# 创建一个空的DataFrame列表来存储所有读取的DataFrame  
df_list = []  
  
# 循环读取所有CSV文件  
for year in range(8, 12):  # 从2008年到2011年  
    filename = fr'C:\Users\zjl15\PycharmProjects\pythonProject1\shixundata\2\acc_{year:02d}_final_mini.csv'  # 注意fr的使用  
    df = pd.read_csv(filename)  # 读取CSV文件  
    df_list.append(df)  # 将读取的DataFrame添加到列表中  
  
# 使用concat函数进行纵向合并（堆叠）  
combined_df = pd.concat(df_list, ignore_index=True)  
  
# 显示合并后的DataFrame  
print(combined_df)  
  
#%%
# #如果你需要将合并后的DataFrame保存到新的CSV文件  
# combined_df.to_csv('combined_final_mini.csv')
#%%
# 将 TRADE_DATE 列转换为 datetime 类型
combined_df['TRADE_DATE'] = pd.to_datetime(combined_df['TRADE_DATE'], format='%Y-%m-%d-%H.%M.%S.%f')

# 使用 apply 函数提取日期部分并创建新的列
combined_df['DATE_ONLY'] = combined_df['TRADE_DATE'].apply(lambda x: x.date())

# 打印结果
print(combined_df)

#%%
df=combined_df.drop(columns=['TRADE_DATE'])
print(df)
#%%
# 使用 groupby 和 transform 计算每个站点每天的交易数量并将结果作为新列加入原数据
df['COUNT'] = df.groupby(['DATE_ONLY', 'TRADE_ADDRESS'])['TRADE_TYPE'].transform('count')

# 打印结果
print(df)
#%%
df=df.drop(columns=['TRADE_TYPE'])
print(df)
#%% md

#%%
# 按照 TRADE_ADDRESS 和 DATE_ONLY 进行排序
df_sorted = df.sort_values(by=['TRADE_ADDRESS', 'DATE_ONLY'])

# 去重处理
df_unique = df_sorted.drop_duplicates(subset=['TRADE_ADDRESS', 'DATE_ONLY'])

# 打印结果
print(df_unique)

#%%
# 根据日期和站点对数据进行分组，并计算客流量总和
grouped_data = df_unique.groupby(['DATE_ONLY', 'TRADE_ADDRESS'])['COUNT'].sum().unstack()

# 绘制折线图
grouped_data.plot(figsize=(12, 8))
plt.title('Daily Traffic for Each Station')
plt.xlabel('Date')
plt.ylabel('Traffic Count')
plt.legend(title='Station', bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

#%%
# 从原始数据中筛选出TRADE_ADDRESS=121的行
data_121 = df_unique[df_unique['TRADE_ADDRESS'] == 121]
data = data_121
print(data_121)
data_121 = data_121.drop(columns=['TRADE_ADDRESS','DATE_ONLY'])
print(data_121)
#%%
# 划分训练集和测试集
train_set = data_121.iloc[:115]
test_set = data_121.iloc[115:]

print("训练集:")
print(train_set)
print("\n测试集:")
print(test_set)

#%%
# 绘制时序图
plt.figure(figsize=(12, 6))
plt.plot(train_set['COUNT'])
plt.title('Traffic Count at Station 121 (Training Set)')
plt.xlabel('Date')
plt.ylabel('Traffic Count')
plt.show()

# 绘制自相关图
fig, ax = plt.subplots(2, 1, figsize=(12, 8))
sm.graphics.tsa.plot_acf(train_set['COUNT'], lags=30, ax=ax[0])
sm.graphics.tsa.plot_pacf(train_set['COUNT'], lags=30, ax=ax[1])
plt.show()

# ADF单位根检测
result = adfuller(train_set['COUNT'])
print('ADF Statistic:', result[0])
print('p-value:', result[1])
print('Critical Values:', result[4])

#%%
# 对训练集进行2阶差分处理
train_set_diff = train_set['COUNT'].diff().diff().dropna()

# 绘制2阶差分后的时序图
plt.figure(figsize=(12, 6))
plt.plot(train_set_diff)
plt.title('Differenced Traffic Count at Station 121 (Training Set)')
plt.xlabel('Date')
plt.ylabel('Differenced Traffic Count')
plt.show()

# 绘制2阶差分后的自相关图
fig, ax = plt.subplots(2, 1, figsize=(12, 8))
sm.graphics.tsa.plot_acf(train_set_diff, lags=30, ax=ax[0])
sm.graphics.tsa.plot_pacf(train_set_diff, lags=30, ax=ax[1])
plt.show()

# ADF单位根检测（2阶差分后）
result_diff = adfuller(train_set_diff)
print('ADF Statistic (after differencing):', result_diff[0])
print('p-value (after differencing):', result_diff[1])
print('Critical Values (after differencing):', result_diff[4])

#%%
# 计算残差序列
residuals = train_set_diff # 这里填入差分处理后的训练集数据的残差序列
# 进行Ljung-Box检验
lb_test_stat, lb_p_value = acorr_ljungbox(residuals, lags=1)
# 输出检验结果
print("Ljung-Box test statistic:", lb_test_stat)
print("P-values:", lb_p_value)
#%%
# 定义 ARIMA 模型的参数 (p, d, q)
p = 6  # 自回归阶数
d = 2  # 差分阶数
q = 6  # 移动平均阶数

# 训练 ARIMA 模型
model = ARIMA(train_set, order=(p, d, q))
model_fit = model.fit()

print(model_fit.summary())

# 使用训练好的模型进行预测
forecast = model_fit.predict(typ='levels')
print(forecast)
# 绘制训练集真实值和预测值对比图
plt.figure(figsize=(12, 6))
plt.plot(train_set, label='Actual')
plt.plot(forecast, label='Predicted', color='red')
plt.title('Comparison of Actual and Predicted Values')
plt.xlabel('Time')
plt.ylabel('Value')
plt.legend()
plt.show()

#%% md

#%%
y_arima = model_fit.forecast(7)
print('预测未来'+str(len(y_arima))+'天的客流量:\n',y_arima)
y_arima.index = test_set.index
# # # 创建一个新的DataFrame来保存预测值，索引与测试集相同  
# test_set = pd.read_csv(r'test.csv')
# print(test_set)
# test_set.to_csv('test_set.xlsx', index=False)
# y_arima = pd.Series(y_arima, index=test_set.index)  
test_set['p'] = y_arima
plt.figure()  
plt.plot(test_set, label='Actual')  
# plt.plot(y_arima, label='Predicted', color='red')  

plt.legend(['y_true','y_arima'])  
plt.show()

#%%
# （2）删掉TRADE_ADDRESS列
df = data.drop(columns=['TRADE_ADDRESS'])
print(df)
#%%
# 将 DATE_ONLY 列转换为日期时间类型
df['DATE_ONLY'] = pd.to_datetime(df['DATE_ONLY'])

# 使用 apply 函数计算星期几，并赋值给新列 DAY_OF_WEEK
df['DAY_OF_WEEK'] = df['DATE_ONLY'].dt.day_name()

# 使用 dt.dayofweek + 1 将星期几转换为数字表示，并赋值给新列 DAY_OF_WEEK_NUM
df['DAY_OF_WEEK_NUM'] = df['DATE_ONLY'].dt.dayofweek + 1

# 显示结果
print(df)
#%%
df2 = df.drop(columns=['DAY_OF_WEEK'])

print(df2)
#%%
# 将 TRADE_DATE 列转换为日期时间类型
df2['DATE_ONLY'] = pd.to_datetime(df2['DATE_ONLY'])

# 使用 apply 函数结合 dt.month 提取出月份信息，并赋值给新列 MONTH
df2['MONTH'] = df2['DATE_ONLY'].apply(lambda x: x.month)
# 使用 apply 函数结合 dt.day 提取出日信息，并赋值给新列 DAY
df2['DAY'] = df2['DATE_ONLY'].apply(lambda x: x.day)
# 显示结果
print(df2)
#%%
df2 = df2.drop(columns=['DATE_ONLY'])
print(df2)
#%%
# 自定义函数来判断是否为周末
def is_weekend(day_of_week):
    return 1 if day_of_week >= 6 else 0

# 使用 apply 函数调用自定义函数对 DAY_OF_WEEK_NUM 列进行计算，生成 WEEKEND 列
df2['WEEKEND'] = df2['DAY_OF_WEEK_NUM'].apply(is_weekend)

print(df2)
#%%
df2 = df2[['MONTH', 'DAY', 'DAY_OF_WEEK_NUM', 'WEEKEND', 'COUNT']]
print(df2)
#%%
# 绘制多变量联合分布图
sns.pairplot(df2)
plt.show()
#%%
from sklearn.model_selection import train_test_split
# 提取特征 x
X = df2[['MONTH', 'DAY', 'DAY_OF_WEEK_NUM', 'WEEKEND']]  # 选择需要作为特征的列

# 提取目标值 y
y = df2['COUNT']  # 选择目标列

# 划分训练集和测试集
# 划分训练集和测试集
X_train = X[:115]  # 前115条作为训练集
X_test = X[115:]   # 后7条作为测试集
y_train = y[:115]  # 对应的目标值作为训练集
y_test = y[115:]   # 对应的目标值作为测试集
print(X_train)
print(y_train)
print(X_test)
print(y_test)
#%%
#构建模型
model = tf.keras.Sequential([
    tf.keras.layers.Dense(64,activation='relu',input_shape=[len(X_train.keys())]),
    tf.keras.layers.Dropout(0.3),
    tf.keras.layers.Dense(64,activation='relu'),
    tf.keras.layers.Dropout(0.3),
    tf.keras.layers.Dense(64,activation='relu'),
    tf.keras.layers.Dropout(0.3),
    tf.keras.layers.Dense(1)
])

# 编译模型
model.compile(optimizer='adam', loss='mse', metrics=['mae','mse'])

model.summary()  

# 训练模型并保存历史数据
history = model.fit(X_train, y_train, epochs=3000, validation_split=0.3)

#%%
#model.save('roch_classification2_cnn.keras')
#%%
# 获取训练过程中的指标数值
mae = history.history['mae']
mse = history.history['loss']
val_mae = history.history['val_mae']
val_mse = history.history['val_mse']

plt.plot(mae, label='MAE')
plt.plot(val_mae, label='Validation MAE')
plt.title('MAE vs')
plt.legend()
plt.show()
plt.plot(mse, label='MSE')
plt.plot(val_mse, label='Validation MSE')
plt.title('MSE')
plt.legend()
plt.show()

# history.plot()
# plt.plot(history)
# plt.show()
print('MAE:', mae)
print('MSE:', mse)
#%%
y_mlp = model.predict(X_test)
#%%
y_test = pd.DataFrame(y_test)
y_mlp = pd.DataFrame(y_mlp)
y_mlp.index = y_arima.index
y_test['P1']=y_arima
y_test['P2']=y_mlp
print(y_test)
plt.figure()  
plt.plot(y_test, label='Actual')  
plt.legend(['y_true','y_arima','y_mlp'])  
plt.show()
#%%
# 计算MAE,MSE,RMSE,MAPE
from sklearn import metrics
import math
MAE = metrics.mean_absolute_error(X_test.iloc[:,0],X_test.iloc[:,1])
MSE = metrics.mean_squared_error(X_test.iloc[:,0],X_test.iloc[:,1])
RMSE = math.sqrt(MSE)
MAPE = metrics.mean_absolute_percentage_error(X_test.iloc[:,0],X_test.iloc[:,1])
print('MAE:',MAE)
print('MSE:',MSE)
print('RMSE:',RMSE)
print('MAPE:',MAPE)
#%%
from sklearn.ensemble import RandomForestRegressor
from sklearn.datasets import make_regression
from sklearn.model_selection import train_test_split
from sklearn import metrics
from sklearn.metrics import mean_squared_error
# 初始化随机森林回归器
rf_regressor = RandomForestRegressor(n_estimators=2500,max_depth=100,random_state=42)
# 训练随机森林模型
rf_regressor.fit(X_train, y_train)
# 使用训练好的模型进行预测
y_regressor_pred = rf_regressor.predict(X_test)
print(y_regressor_pred)
y_regressor_pred = pd.DataFrame(y_regressor_pred)
print(y_regressor_pred)
y_regressor_pred.index = y_mlp.index
y_test['P3'] = y_regressor_pred
print(y_test)
plt.figure()  
plt.plot(y_test, label='Actual')  
plt.legend(['y_true','y_arima','y_mlp','y_regressor_pred'])  
plt.show()
# 计算均方误差
mse = mean_squared_error(y_test['COUNT'], y_test['P3'])
print("Mean Squared Error:", mse)
#%%
from xgboost import XGBRegressor
from sklearn.metrics import mean_squared_error
# 1. 准备数据集
# 假设数据集已准备好，包含特征和目标变量
# 4. 定义模型
model = XGBRegressor(
    n_estimators=3000,  # 树的数量
    max_depth=5,       # 树的最大深度
    learning_rate=0.3, # 学习率
    random_state=42
)
# 5. 训练模型
model.fit(X_train, y_train)
# 6. 模型评估
y_XGBRegressor_pred = model.predict(X_test)
y_XGBRegressor_pred = pd.DataFrame(y_XGBRegressor_pred)
y_XGBRegressor_pred.index = y_mlp.index
y_test['P4'] = y_XGBRegressor_pred
print(y_test)
plt.figure()  
plt.plot(y_test, label='Actual')  
plt.legend(['y_true','y_arima','y_mlp','y_regressor_pred','y_XGBRegressor_pred'])  
plt.show()
# 计算均方误差
mse = mean_squared_error(y_test['COUNT'], y_test['P3'])
print("Mean Squared Error:", mse)
mse = mean_squared_error(y_test['COUNT'], y_XGBRegressor_pred)
print("Mean Squared Error:", mse)
#%%
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.svm import SVR
from sklearn.metrics import mean_squared_error

# 4. 定义模型
model = SVR(kernel='rbf', C=1.0, epsilon=0.1)

# 5. 训练模型
model.fit(X_train, y_train)

# 6. 模型评估
y_SVR_pred = model.predict(X_test)
y_SVR_pred = pd.DataFrame(y_SVR_pred)
y_SVR_pred.index = y_mlp.index
y_test['P5'] = y_SVR_pred
print(y_test)
plt.figure()  
plt.plot(y_test, label='Actual')  
plt.legend(['y_true','y_arima','y_mlp','y_regressor_pred','y_XGBRegressor_pred','y_SVR_pred'])  
plt.show()
# 计算均方误差
mse = mean_squared_error(y_test['COUNT'], y_test['P5'])
print("Mean Squared Error:", mse)
mse = mean_squared_error(y_test['COUNT'], y_SVR_pred)
print("Mean Squared Error:", mse)
#%%