机器学习之心【Python机器学习/深度学习程序全家桶】

#!/usr/bin/env python# coding: utf-8
# In[1]:

# 导入环境中的相关包import itertools  # 导入迭代工具模块import numpy as np  # 导入NumPy库import pandas as pd  # 导入Pandas库import matplotlib.pyplot as plt  # 导入Matplotlib库from matplotlib.ticker import MultipleLocator  # 导入刻度定位器import warnings  # 导入警告模块from statsmodels.graphics.tsaplots import plot_acf, plot_pacf  # 导入自相关图和偏自相关图绘制函数from statsmodels.stats.diagnostic import acorr_ljungbox  # 导入Ljung-Box检验函数from statsmodels.tsa.statespace.sarimax import SARIMAX  # 导入SARIMAX模型from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error  # 导入评估指标函数from statsmodels.tsa.stattools import adfuller  # 导入ADF检验函数import math  # 导入数学函数库#更多模型咸鱼搜索机器学习之心，支持模型定制 import seaborn as sns  # 导入Seaborn库import statsmodels.api as sm  # 导入Statsmodels库import tensorflow as tf  # 导入TensorFlow库from pmdarima import auto_arima  # 导入自动ARIMA模型选择函数from statsmodels.tsa.seasonal import seasonal_decompose  #导入分解函数from scipy.signal import find_peaks  #用于查找信号中峰值（波峰）位置的函数

# In[2]:

# 显示中文plt.rcParams['font.sans-serif'] = ['SimHei']  # 设置字体为SimHeiplt.rcParams['figure.figsize'] = (10.0, 8.0)  # 设置默认图像大小为10x8英寸plt.rcParams['image.interpolation'] = 'nearest'  # 设置图像插值方式为最近邻插值plt.rcParams['image.cmap'] = 'gray'  # 设置图像颜色映射为灰度图#更多模型咸鱼搜索机器学习之心，支持模型定制 

# In[3]:

# 忽略警告warnings.filterwarnings('ignore')

# In[4]:

#更多模型咸鱼搜索机器学习之心，支持模型定制 

# In[5]:

df = pd.read_csv("A.csv", header=0, parse_dates=[0], index_col=0, usecols=[0, 1]) #读取数据df.head()#显示数据前5行

# In[6]:

# 创建一个分辨率为150点每英寸（dpi），尺寸为15x3英寸的图形和轴对象#更多模型咸鱼搜索机器学习之心，支持模型定制 fig, ax = plt.subplots(dpi=150, figsize=(15, 3))# 使用线宽为1，在'ax'轴对象上绘制'dataframe'中的数据df.plot(ax=ax, linewidth=1)# 设置图形的标题为'A'ax.set_title('A')# 设置y轴的标签为'A'ax.set_ylabel('A')# 设置x轴的标签为'time/day'ax.set_xlabel('time/day')# 设置字体样式为'Times New Roman'，字体大小为10font = {'serif': 'Times New Roman', 'size': 10}plt.rc('font', **font)# 显示图形plt.show()

# In[7]:

data = df['A'].values  # 计算自相关系数correlation = np.correlate(data, data, mode='full')

# In[8]:

# 绘制自相关图plt.plot(correlation)#更多模型咸鱼搜索机器学习之心，支持模型定制 plt.title('自相关图')plt.xlabel('滞后')plt.ylabel('相关性')plt.show()

# In[9]:

# 找到波峰的位置并计算波峰之间的距离peaks, _ = find_peaks(correlation)peak_distances = np.diff(peaks)estimated_period = np.mean(peak_distances)
print("估计周期为:", estimated_period)

# In[10]:

# 进行季节性分解result = seasonal_decompose(df['A'], model='multiplicative', period=int(estimated_period))  
# 绘制分解图plt.figure(figsize=(12, 8))plt.subplot(411)plt.plot(df['A'], label='Original')plt.legend(loc='upper left')#更多模型咸鱼搜索机器学习之心，支持模型定制 plt.subplot(412)plt.plot(result.trend, label='Trend')plt.legend(loc='upper left')
plt.subplot(413)plt.plot(result.seasonal, label='Seasonal')plt.legend(loc='upper left')#更多模型咸鱼搜索机器学习之心，支持模型定制 plt.subplot(414)plt.plot(result.resid, label='Residual')plt.legend(loc='upper left')
plt.tight_layout()plt.show()

# In[11]:

#定义稳定性检验函数def adf_val(ts, ts_title):    # 计算ADF检验值    adf, pvalue, usedlag, nobs, critical_values, icbest = adfuller(ts)    # 定义指标的名称和对应的值    name = ['adf', 'pvalue', 'usedlag', 'nobs', 'critical_values', 'icbest']    values = [adf, pvalue, usedlag, nobs, critical_values, icbest]    # 打印各项指标的名称和对应的值    print(list(zip(name, values)))    # 返回ADF值、p值和临界值    return adf, pvalue, critical_values    # 返回adf值、adf的p值、三种状态的检验值    #更多模型咸鱼搜索机器学习之心，支持模型定制 

# In[12]:

#白噪声检验也称为纯随机性检验，当数据是纯随机数据时，再对数据进行分析就没有任何意义了，所以拿到数据后最好对数据进行一个纯随机性检验。def acorr_val(ts):    '''    # 白噪声（随机性）检验    ts: 时间序列数据，Series类型    返回白噪声检验的P值    '''    #更多模型咸鱼搜索机器学习之心，支持模型定制     lb_test_results = acorr_ljungbox(ts, lags=2)  # 白噪声检验结果
    return lb_test_results # 返回白噪声检验的统计量值和P值

# In[13]:

def tsplot(y, lags=None, figsize=(14, 8)):    fig = plt.figure(figsize=figsize)  # 创建绘图窗口，设置尺寸为figsize    layout = (2, 2)  # 子图布局为2行2列    ts_ax = plt.subplot2grid(layout, (0, 0))  # 创建时间序列图的子图对象    hist_ax = plt.subplot2grid(layout, (0, 1))  # 创建直方图的子图对象    acf_ax = plt.subplot2grid(layout, (1, 0))  # 创建自相关图的子图对象    pacf_ax = plt.subplot2grid(layout, (1, 1))  # 创建偏自相关图的子图对象
    y.plot(ax=ts_ax)  # 绘制时间序列图    ts_ax.set_title('A Given Training Series')  # 设置时间序列图标题
    y.plot(ax=hist_ax, kind='hist', bins=25)  # 绘制直方图    hist_ax.set_title('Histogram')  # 设置直方图标题
    plot_acf(y, lags=lags, ax=acf_ax)  # 绘制自相关图    plot_pacf(y, lags=lags, ax=pacf_ax)  # 绘制偏自相关图
    [ax.set_xlim(0) for ax in [acf_ax, pacf_ax]]  # 设置自相关图和偏自相关图x轴刻度起点为0
    sns.despine()  # 移除图形顶部和右侧的边框
    fig.tight_layout()  # 调整子图布局和间距    fig.show()  # 显示图形
    return ts_ax, acf_ax, pacf_ax  # 返回坐标轴对象#更多模型咸鱼搜索机器学习之心，支持模型定制 

# In[14]:

ts_data = df.astype('float32')  # 将DataFrame df的数据类型转换为'float32'类型，以确保时间序列数据为浮点数类型

# In[15]:

#adf结果为-10.4， 小于三个level的统计值。pvalue也是接近于0 的，所以是平稳的#更多模型咸鱼搜索机器学习之心，支持模型定制 adf, pvalue1, critical_values = adf_val(ts_data, 'raw time series')print('adf',adf)print('pvalue1',pvalue1)print('critical_values',critical_values)#若p值远小于0.01，认为该时间序列是平稳的aco=acorr_val(ts_data)print('aco',aco)

# In[16]:

diff1 = ts_data.diff()#对ts_data数据框进行一阶差分操作，即每个时间点的值减去前一个时间点的值，得到的差分结果赋值给变量diff1。#画出一阶差分后的数据#更多模型咸鱼搜索机器学习之心，支持模型定制 plt.figure(figsize=(12,3))plt.plot(diff1, 'blue', label='diff1')plt.title('diff1')plt.xlabel('time')plt.ylabel('price')plt.legend()plt.show()

# In[17]:

'''将diff1中的缺失值（NaN）用0进行填充，得到填充后的差分结果。然后，将差分结果中的正无穷和负无穷值（Inf）用0进行替换，得到最终处理后的差分结果。这两行代码的目的是处理差分结果中可能出现的异常值。'''diff1[np.isnan(diff1)] = 0diff1[np.isinf(diff1)] = 0

# In[18]:

#对差分后的数据再次进行检验adf, pvalue1, critical_values = adf_val(diff1, 'raw time series')print('adf',adf)print('pvalue1',pvalue1)print('critical_values',critical_values)aco=acorr_val(ts_data)print('aco',aco)

# In[19]:

diff2 = ts_data.diff().diff()#对ts_data数据框进行二阶差分操作。#更多模型咸鱼搜索机器学习之心，支持模型定制 plt.figure(figsize=(12,3))plt.plot(diff2, 'r', label='diff2')plt.title('diff2')plt.xlabel('time')plt.ylabel('price')plt.legend()plt.show()

# In[20]:

'''将diff2中的缺失值（NaN）用0进行填充，得到填充后的差分结果。然后，将差分结果中的正无穷和负无穷值（Inf）用0进行替换，得到最终处理后的差分结果。这两行代码的目的是处理差分结果中可能出现的异常值。'''diff2[np.isnan(diff2)] = 0diff2[np.isinf(diff2)] = 0

# In[21]:

#对二次差分结果进行检验adf, pvalue1, critical_values = adf_val(diff2, 'raw time series')print('adf',adf)print('pvalue1',pvalue1)print('critical_values',critical_values)aco=acorr_val(ts_data)print('aco',aco)#更多模型咸鱼搜索机器学习之心，支持模型定制 

# In[22]:

##自相关和偏自相关tsplot(ts_data, lags=20)#更多模型咸鱼搜索机器学习之心，支持模型定制 

# In[23]:

train_data, test_data = df[0:int(len(df)*0.8)], df[int(len(df)*0.8):]  # 划分训练集和测试集，按照80:20比例

# In[24]:

#取划分的数据train_ar = train_data.values  # 将训练集数据转换为 NumPy 数组test_ar = test_data.values  # 将测试集数据转换为 NumPy 数组

# In[ ]:




# '''# 使用 auto_arima 函数对训练集数据 train_data 进行自动ARIMA模型拟合。# 设置 seasonal=True 表示考虑季节性，m=12 表示季节性周期为12个时间步长。# max_p=7、max_d=2、max_q=7 分别表示自回归、差分和移动平均部分的最大阶数。# max_P=4、max_D=4、max_Q=4 分别表示季节性自回归、差分和移动平均部分的最大阶数。# 使用 .summary() 方法输出ARIMA模型的摘要信息。# '''# auto_arima(train_data, seasonal=True, m=12,max_p=7, max_d=2,max_q=7, max_P=4, max_D=4,max_Q=4).summary()
# In[ ]:




# In[25]:

def best_sarima_model(train_data, p, q, P, Q, d=1, D=1, s=12):    best_model_aic = np.Inf  # 初始化最佳模型的 AIC 值为正无穷大    best_model_bic = np.Inf  # 初始化最佳模型的 BIC 值为正无穷大    best_model_hqic = np.Inf  # 初始化最佳模型的 HQIC 值为正无穷大    best_model_order = (0, 0, 0)  # 初始化最佳模型的阶数    models = []  # 用于存储所有拟合的模型    #更多模型咸鱼搜索机器学习之心，支持模型定制     # 遍历所有可能的阶数组合    for p_ in p:        for q_ in q:            for P_ in P:                for Q_ in Q:#更多模型咸鱼搜索机器学习之心，支持模型定制                     try:                        no_of_lower_metrics = 0  # 记录优于当前最佳模型指标的个数                        model = SARIMAX(endog=train_data, order=(p_, d, q_), seasonal_order=(P_, D, Q_, s),                                        enforce_invertibility=False).fit()  # 拟合 SARIMA 模型                        models.append(model)  # 存储拟合的模型
                        # 计算模型指标并与当前最佳模型进行比较                        if model.aic <= best_model_aic:                            no_of_lower_metrics += 1                        if model.bic <= best_model_bic:                            no_of_lower_metrics += 1                        if model.hqic <= best_model_hqic:#更多模型咸鱼搜索机器学习之心，支持模型定制                             no_of_lower_metrics += 1
                        # 如果有至少两个指标优于当前最佳模型，则更新最佳模型的指标和阶数                        if no_of_lower_metrics >= 2:                            best_model_aic = np.round(model.aic, 0)                            best_model_bic = np.round(model.bic, 0)                            best_model_hqic = np.round(model.hqic, 0)                            best_model_order = (p_, d, q_, P_, D, Q_, s)                            current_best_model = model                            models.append(model)#更多模型咸鱼搜索机器学习之心，支持模型定制                             print("Best model: SARIMA" +  str(best_model_order) +                                   " AIC:{} BIC:{} HQIC:{}".format(best_model_aic, best_model_bic, best_model_hqic) +                                  " resid:{}".format(np.round(np.exp(current_best_model.resid).mean(), 3)))
                    except:                        pass
    print('\n')    print(current_best_model.summary())  # 输出最佳模型的详细信息    return current_best_model, models
# 调用 best_sarima_model 函数，传入训练数据 train_ar 和参数候选值# 将返回的最佳模型赋值给 best_model，将所有拟合的模型列表赋值给 modelsbest_model, models = best_sarima_model(train_data=train_ar, p=range(3), q=range(3), P=range(3), Q=range(3))

# In[ ]:




# In[ ]:

# 定义 p、d、q 的取值范围p = range(0, 3)d = range(0, 1)q = range(0, 3)#更多模型咸鱼搜索机器学习之心，支持模型定制 #更多模型咸鱼搜索机器学习之心，支持模型定制 # 生成 pdq 参数组合pdq = list(itertools.product(p, d, q))# 生成 seasonal_pdq 参数组合seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]#更多模型咸鱼搜索机器学习之心，支持模型定制 min_aic = 999999999
# 遍历 pdq 参数组合和 seasonal_pdq 参数组合for param in pdq:    for param_seasonal in seasonal_pdq:        try:            # 创建 SARIMAX 模型并拟合训练数据            mod = sm.tsa.statespace.SARIMAX(train_ar,                                            order=param,                                            seasonal_order=param_seasonal,                                            enforce_stationarity=False,                                            enforce_invertibility=False)            results = mod.fit()
            # 打印模型的 AIC 值            print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
            # 更新最小 AIC 值和对应的模型            if results.aic < min_aic:                min_aic = results.aic                min_aic_model = results
        except:            continue#更多模型咸鱼搜索机器学习之心，支持模型定制 # 打印具有最小 AIC 值的模型的摘要信息print(min_aic_model.summary())

# In[ ]:




# In[26]:

# 构建训练数据history = [x for x in train_ar]  # 使用训练数据初始化历史数据列表predictions = list()  # 初始化预测结果列表#更多模型咸鱼搜索机器学习之心，支持模型定制 # 训练ARIMA模型并进行预测for t in range(len(test_ar)):  # 遍历测试数据集的每个时间步    model = sm.tsa.SARIMAX(history, order=(2, 1, 1), seasonal_order=(2, 1, 2, 12), enforce_invertibility=False)  # 构建SARIMAX模型    model_fit = model.fit()  # 拟合模型    output = model_fit.forecast()  # 模型预测    yhat = output[0]  # 获取预测结果    #更多模型咸搜索机器学习之心，支持模型定制     predictions.append(yhat)  # 将预测结果添加到预测结果列表中    obs = test_ar[t]  # 获取当前观测值    history.append(obs)  # 将当前观测值添加到历史数据列表中    print('predicted=%f, expected=%f' % (yhat, obs))  # 打印预测值和真实值

# In[27]:

testScore = math.sqrt(mean_squared_error(test_ar, predictions))  # 计算均方根误差（RMSE）print('RMSE %.3f' % testScore)  # 打印RMSE值#更多模型咸鱼搜索机器学习之心，支持模型定制 testScore = r2_score(test_ar, predictions)  # 计算R平方值（R2）print('R2 %.3f' % testScore)  # 打印R2值#更多模型咸鱼搜索机器学习之心，支持模型定制 testScore = mean_absolute_error(test_ar, predictions)  # 计算平均绝对误差（MAE）print('MAE %.3f' % testScore)  # 打印MAE值

# In[28]:

#只显示预测部分，不显示训练数据部分plt.figure(figsize=(12, 7))  # 创建绘图窗口，设置尺寸
plt.plot(test_data.index, predictions, color='b', marker='o', linestyle='dashed', label='Predicted')  # 绘制预测结果曲线plt.plot(test_data.index, test_data, color='red', label='Actual')  # 绘制真实值曲线
plt.title('SARIMA')  # 设置图标题plt.xlabel('time')  # 设置x轴标签plt.ylabel('A')  # 设置y轴标签plt.legend()  # 显示图例
plt.show()  # 显示图形

# In[ ]:




# In[ ]: