基于时间序列的温度预测

一、选题背景
   人们的一切活动，其目的无不在认识世界和改造世界，时间序列分析不仅可以从数量上揭示某一现象的发展变化规律或从动态
的角度刻划某一现象之间与其他现象之间的内在数量关系及其变化规律性，达到认识客观世界之目的。而且运用时序模型还可以
预测和控制现象的未来行为，修正或重新设计系统以达到利用和改造客观世界之目的。近几年来，时间序列分析引起了国内外学
者及科研和管理人员的极大兴趣，特别是随着计算机的普及和软件的开发应用，对于只具有一般数学知识的学者和广大的工程
技术及管理人员学习和掌握时间序列分析方法，并用以分析、探索社会经济现象的动态结构和发展变动规律，进行对未来状态进
行预测控制，提供了实现可能性，且在诸多应用领域已取得了可喜成果。
   温度对一个地区的农业，工业，生活都具有很重要的意义，温度随时间的变化而不同，其所形成的序列可以看成是时间序列。

二、机器学习案例设计方案
1.数据集来源
kaggle耶拿天气数据集
2009年1月10日至2016年12月31日，每10分钟记录一次，具有14个特征，例如温度，压力，湿度等。


2.机器学习框架描述
Apache Singa 是一个用于在大型数据集上训练深度学习的通用分布式深度学习平台，它是基于分层抽象的简单开发模型设计
的。它还支持各种当前流行的深度学习模型，有前馈模型（卷积神经网络，CNN），能量模型（受限玻尔兹曼机，RBM和循环神
经网络，RNN），还为用户提供了许多内嵌层。

3.遇到的困难

数据预处理
构建模型
训练集，验证集



import pandas as pd
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow import keras

csv_path = "jena_climate_2009_2016.csv"
df = pd.read_csv(csv_path)

#标题
titles = [
    "Pressure",
    "Temperature",
    "Temperature in Kelvin",
    "Temperature (dew point)",
    "Relative Humidity",
    "Saturation vapor pressure",
    "Vapor pressure",
    "Vapor pressure deficit",
    "Specific humidity",
    "Water vapor concentration",
    "Airtight",
    "Wind speed",
    "Maximum wind speed",
    "Wind direction in degrees",
]

#列名
feature_keys = [
    "p (mbar)",
    "T (degC)",
    "Tpot (K)",
    "Tdew (degC)",
    "rh (%)",
    "VPmax (mbar)",
    "VPact (mbar)",
    "VPdef (mbar)",
    "sh (g/kg)",
    "H2OC (mmol/mol)",
    "rho (g/m**3)",
    "wv (m/s)",
    "max. wv (m/s)",
    "wd (deg)",
]

colors = [
    "blue",
    "orange",
    "green",
    "red",
    "purple",
    "brown",
    "pink",
    "gray",
    "olive",
    "cyan",
]
date_time_key = "Date Time"

df = pd.read_csv('D:\Download\jena_climate_2009_2016.csv')
df.head()

import os

fname = 'jena_climate_2009_2016.csv'

# 读取数据
f = open(fname)
data = f.read()
f.close()

# 分割数据成每一行，保存
lines = data.split('\n')
# 找到数据头
header = lines[0].split(',')
# 剔除数据头
lines = lines[1:]

print(len(lines))
print(header)

# 将数据转成Numpy数组

import numpy as np

float_data = np.zeros((len(lines),len(header)-1))

for i,line in enumerate(lines):
    values = [float(x) for x in line.split(',')[1:]]
    float_data[i,:] = values
# 绘制温度时间序列

from matplotlib import pyplot as plt

temp = float_data[:,1]
plt.plot(range(len(temp)),temp)

plt.figure() 

plt.plot(range(1440), temp[:1440])

def show_raw_visualization(data):
    time_data = data[date_time_key]
    # 7*2的总图形结构
    fig, axes = plt.subplots(
        nrows=7, ncols=2, figsize=(15, 20), dpi=80, facecolor="w", edgecolor="k"
    )
    for i in range(len(feature_keys)):
        key = feature_keys[i]
        c = colors[i % (len(colors))] # 选择一个颜色，循环选择
        t_data = data[key] # 获取对应字段数据
        t_data.index = time_data # 设置索引为时间
        t_data.head()
        ax = t_data.plot(
            ax=axes[i // 2, i % 2],
            color=c,
            title="{} - {}".format(titles[i], key),
            rot=25,
        )
        ax.legend([titles[i]])
    plt.tight_layout()   # 这个的作用是自适应图片大小

#原始数据可视化
show_raw_visualization(df)

def show_heatmap(data):
    plt.matshow(data.corr())
    plt.xticks(range(data.shape[1]), data.columns, fontsize=14, rotation=90)
    plt.gca().xaxis.tick_bottom()
    plt.yticks(range(data.shape[1]), data.columns, fontsize=14)

    cb = plt.colorbar()
    cb.ax.tick_params(labelsize=14)
    plt.title("Feature Correlation Heatmap", fontsize=14)
    plt.show()
    
# 展示字段间的相关性
show_heatmap(df)

#数据预处理

split_fraction = 0.715 # 训练和测试数据集划分比例
train_split = int(split_fraction * int(df.shape[0])) # 计算出划分的索引处
step = 6 # 数据的选取步频

past = 720 # 使用前720个数据时间进行预测
future = 72 # 预测72个数据时点后的数据
learning_rate = 0.001 # 学习速率，用于建模的优化器参数
batch_size = 256 # 训练的数据批此规模
epochs = 25 # 训练25次

def normalize(data, train_split):
    data_mean = data[:train_split].mean(axis=0)
    data_std = data[:train_split].std(axis=0)
    return (data - data_mean) / data_std

#相对湿度等是多余的

print(
    "The selected parameters are:",
    ", ".join([titles[i] for i in [0, 1, 5, 7, 8, 10, 11]]),
)

# 选择更合适的字段进行建模
selected_features = [feature_keys[i] for i in [0, 1, 5, 7, 8, 10, 11]]
features = df[selected_features]
features.index = df[date_time_key]
features.head()

# 所有数据进行归一化
features = normalize(features.values, train_split)
features = pd.DataFrame(features)
features.head()

# 划分训练数据和测试数据
train_data = features.loc[0 : train_split - 1] # 训练数据
val_data = features.loc[train_split:]  # 测试数据

print(train_data.head(2))
print(val_data.head(2))


start = past + future
end = start + train_split

x_train = train_data[[i for i in range(7)]].values
y_train = features.iloc[start:end][[1]]

sequence_length = int(past / step)

# 构造训练数据
dataset_train = keras.preprocessing.timeseries_dataset_from_array(
    x_train,
    y_train,
    sequence_length=sequence_length, # 一批采集120次
    sampling_rate=step,  # 6步采集一次，即 60分钟采集一次
    batch_size=batch_size,
)

# 构造验证数据

x_end = len(val_data) - past - future  # 因为是用720个监测数据去预测72个窗口数据后的值，所以要倒减一下

label_start = train_split + past + future  # 同理Y值要延后一下

x_val = val_data.iloc[:x_end][[i for i in range(7)]].values
y_val = features.iloc[label_start:][[1]]

dataset_val = keras.preprocessing.timeseries_dataset_from_array(
    x_val,
    y_val,
    sequence_length=sequence_length,
    sampling_rate=step,
    batch_size=batch_size,
)


for batch in dataset_train.take(1):
    inputs, targets = batch

print("Input shape:", inputs.numpy().shape)
print("Target shape:", targets.numpy().shape)

# 创建模型
inputs = keras.layers.Input(shape=(inputs.shape[1], inputs.shape[2]))
lstm_out = keras.layers.LSTM(32)(inputs)
outputs = keras.layers.Dense(1)(lstm_out)

model = keras.Model(inputs=inputs, outputs=outputs)
# optimizer 优化器，加快学习速度，loss 损失计算
model.compile(optimizer=keras.optimizers.Adam(learning_rate=learning_rate), loss="mse")
# 展示模型
model.summary()

# 设置损失函数和训练
path_checkpoint = "model_checkpoint.h5"
es_callback = keras.callbacks.EarlyStopping(monitor="val_loss", min_delta=0, patience=5)

# 使用ModelCheckpoint回调EarlyStopping函数定期保存检查点，并使用该回调函数在验证损失不再改善时中断训练。
modelckpt_callback = keras.callbacks.ModelCheckpoint(
    monitor="val_loss",
    filepath=path_checkpoint,
    verbose=1,
    save_weights_only=True,
    save_best_only=True,
)

history = model.fit(
    dataset_train,
    epochs=epochs,
    validation_data=dataset_val,
    callbacks=[es_callback, modelckpt_callback],
)

# 显示结果
# loss：训练集的差距大小 和 val_loss 测试集的差距大小

def visualize_loss(history, title):
    loss = history.history["loss"]
    val_loss = history.history["val_loss"]
    epochs = range(len(loss))
    plt.figure()
    plt.plot(epochs, loss, "b", label="Training loss")
    plt.plot(epochs, val_loss, "r", label="Validation loss")
    plt.title(title)
    plt.xlabel("Epochs")
    plt.ylabel("Loss")
    plt.legend()
    plt.show()

visualize_loss(history, "Training and Validation Loss")

#模型预测

def show_plot(plot_data, delta, title):
    labels = ["History", "True Future", "Model Prediction"]
    marker = [".-", "rx", "go"]
    time_steps = list(range(-(plot_data[0].shape[0]), 0))
    if delta:
        future = delta
    else:
        future = 0

    plt.title(title)
    for i, val in enumerate(plot_data):
        if i:
             # 绘制预测数据
            plt.plot(future, plot_data[i], marker[i], markersize=10, label=labels[i])
        else:
            # 绘制历史数据
            plt.plot(time_steps, plot_data[i].flatten(), marker[i], label=labels[i])
    plt.legend()
    plt.xlim([time_steps[0], (future + 5) * 2])
    plt.xlabel("Time-Step")
    plt.show()
    return

# 输入5个批此数据，历史数据是120个小时，预测数据是12个小时以后，使用模型预测处结果
for x, y in dataset_val.take(5):
    show_plot(
        [x[0][:, 1].numpy(), y[0].numpy(), model.predict(x)[0]],
        12,
        "Single Step Prediction",
    )