利用Tensorflow2处理时间序列
文章主要参照Andrew Ng与Laurence Moroney的Coursera课程，目的是为记录学习Tensorflow的过程. 文章所有代码内容均使用Python3编写.前言
以下是本篇文章正文内容
一、时间序列
什么是时间序列？
时间序列（或称动态数列）是指将同一统计指标的数值按其发生的时间先后顺序排列而成的数列。时间序列分析的主要目的是根据已有的历史数据对未来进行预测。经济数据中大多数以时间序列的形式给出。根据观察时间的不同，时间序列中的时间可以是年份、季度、月份或其他任何时间形式。
创建简单的时间序列
先导入已有第三方库numpy与matplotlib
import numpy as np
import matplotlib.pyplot as plt
建立两个函数，一个用于绘制序列，另一个用于返回趋势。
def plot_series(time, series):
plt.figure(figsize=(10, 6))
plt.plot(time, series)
plt.xlabel("time")
plt.ylabel("value")
plt.grid(True)
plt.show()
def trend(time, slope=0):
return slope * time
现在，让我们来绘制一个非常简单的时间序列
time = np.arange(4 * 365 + 1)
baseline = 10
series = trend(time, 0.1)
plot_series(time, series)
尽管它是一条直线，但它仍是一个时间序列，这里x轴是时间，y轴是当时函数的值。
时间序列季节性模式
接下来，让我们看一下时间序列的季节性模式。
下面两个函数包含时间序列的季节性模式，然后就是相同模式下的季节性。
def seasonal_pattern(season_time):
return np.where(season_time < 0.4,
np.cos(season_time * 2 * np.pi),
1 / np.exp(3 * season_time))
关于np.where的用法：
np.where(condition,x,y)
满足条件(condition)，输出想，不满足输出y.
在该函数中x即为
c
o
s
2
π
t
cos{2 \pi t}
cos2πt，y即为
e
3
t
e^{3t}
e3t
当然，这里seasonal_pattern只是一个任意的函数，可以人为进行一些调整，得到我们想要的任何形状。
def seasonality(time, period, amplitude=1, phase=0):
season_time = ((time + phase) % period) / period
return amplitude * seasonal_pattern(season_time)
seasonality函数功能是在每个时间段(period)中，重复seasonal_pattern函数的季节性模式.
我们现在将其绘制
baseline = 10
amplitude = 40
series = seasonality(time, period=365, amplitude=amplitude)
plot_series(time, series)
得到下图：
我们通过观察图表，可以看到清晰的波峰和波谷，但除此之外，有较小的有规则的峰值，这可以看成粗略地模拟了时间序列。例如，也许商店的利润在商店关门那天为负值，第二天达到峰值，然后逐渐减小，又在周末达到峰值。
如果这是我们对其添加趋势以至于季节性数据，季节性模式随着时间增加，也许是一个成上升趋势的业务，我们看到相同模式在重复中呈上升趋势。
slope = 0.05
series = baseline + trend(time, slope) + seasonality(time, period=365, amplitude=amplitude)
plot_series(time, series)
这时如果我们再增加一个在时间序列当中常见的因素----噪声，我们来进行进一步观察。建立一个函数添加随机噪声，设置不同的数据，我们就得到比较嘈杂的不同时间序列。
我们还可以通过调整noise函数中的参数noise_level，从而调整曲线的嘈杂程度。
def noise(time, noise_level=1):
return np.random.randn(len(time)) * noise_level
noise_level = 15
noisy_series = series + noise(time, noise_level)
plot_series(time, noisy_series)
noise_level = 40
noisy_series = series + noise(time, noise_level)
plot_series(time, noisy_series)
然而仅通过人眼很难观察出其季节性模式，但计算机能够很好的帮我们发现并解决这个问题。
自相关函数
接下来，将讨论一些自相关的相关内容，有两个不同的自相关函数，接下来将同时绘制这两种图形，以便于看到每种图形的效果。
何为自相关函数？
自相关（英语：Autocorrelation），也叫序列相关， 是一个信号于其自身在不同时间点的互相关。非正式地来说，它就是两次观察之间的相似度对它们之间的时间差的函数。它是找出重复模式（如被噪声掩盖的周期信号），或识别隐含在信号谐波频率中消失的基频的数学工具。它常用于信号处理中，用来分析函数或一系列值，如时域信号。
自相关函数可以理解为序列不同时刻（或不同位置）的相似程度的度量。
下面我们将两个不同自相关函数产生的不同效果分别作图，从而能够清楚地看清其不同的影响。
自相关函数一
def autocorrelation(time, amplitude):
rho1 = 0.5
rho2 = -0.1
ar = np.random.randn(len(time) + 50)
ar[:50] = 100
for step in range(50, len(time) + 50):
ar[step] += rho1 * ar[step - 50]
ar[step] += rho2 * ar[step - 33]
return ar[50:] * amplitude
series = autocorrelation(time, 10)
plot_series(time[:200], series[:200])
series = noise(time)
plot_series(time[:200], series[:200])
series = autocorrelation(time, 10) + trend(time, 2)
plot_series(time[:200], series[:200])
series = autocorrelation(time, 10) + seasonality(time, period=50, amplitude=150) + trend(time, 2)
plot_series(time[:200], series[:200])
series = autocorrelation(time, 10) + seasonality(time, period=50, amplitude=150) + trend(time, 2)
series2 = autocorrelation(time, 5) + seasonality(time, period=50, amplitude=2) + trend(time, -1) + 550
series[200:] = series2[200:]
#series += noise(time, 30)
plot_series(time[:300], series[:300])
自相关函数二
def autocorrelation(time, amplitude):
rho = 0.8
ar = np.random.randn(len(time) + 1)
for step in range(1, len(time) + 1):
ar[step] += rho * ar[step - 1]
return ar[1:] * amplitude
series = autocorrelation(time, 10)
plot_series(time[:200], series[:200])
series = noise(time)
plot_series(time[:200], series[:200])
series = autocorrelation(time, 10) + trend(time, 2)
plot_series(time[:200], series[:200])
series = autocorrelation(time, 10) + seasonality(time, period=50, amplitude=150) + trend(time, 2)
plot_series(time[:200], series[:200])
series = autocorrelation(time, 10) + seasonality(time, period=50, amplitude=150) + trend(time, 2)
series2 = autocorrelation(time, 5) + seasonality(time, period=50, amplitude=2) + trend(time, -1) + 550
series[200:] = series2[200:]
#series += noise(time, 30)
plot_series(time[:300], series[:300])
二、分析与预测
训练集、验证集与测试集
在上面，我们已经了解了所有构成时间序列的不同因素。现在，我们开始研究技术，结合已知，进行对时间序列的预测，我们从时间序列，趋势季节性以及噪声进行分析，从而得到对时间序列的分析结果。
为了衡量预测时间序列模型的效果，我们将时间序列划分为训练期、验证期以及训练期，这些成为固定分期，如果时间序列包含季节性，我们通常希望每个时期都包含整个季节，比如一年两年或者是一个月两个月从而观察时间序列是否具有 年度季节性或者月度季节性等等。
接下来，将在训练期进行训练模型，在验证期进行评估模型，在测试期间测试数据，以查看模型是否表现良好。
如果这样做的话，为了获得更恰当的时间序列模型，我们还可以采取不寻常的方法，即使用测试集进行再次训练，因为测试集是我们所拥有的最接近当前时间的数据，它通常是确定未来数据值的最强依据，如果模型没有使用该数据再训练，则得到的结果很有可能不是最佳结果。
各项指标(Metrics for evaluating performance)
一旦有了模型，我们需对其评估，然后我们需要一个指标计算其效果，因此，我们从计算误差开始。
下面列出了一系列计算误差的式子，包括
- errors:即预测值与验证期间的实际值之差。
- mse:均方误差，即errors平方和的平均值。（使用平方的最主要原因是要消除负值，例如，如果有连个误差值，分别为2和-2，那么使用普通的那种误差，这些误差就会相互抵消，但这显然是错误的，因为我们有两个误差而不是没有，如果对其进行平方，我们就得到两个相等的均方误差）
- rmse:如果希望误差与原来有相同的标度，即不要相差过大，我们就采用均方根误差进行评估
- mae:平均绝对误差，即可以防止负数正数互相抵消，又可以不像mse那样产生较大的误差。
误差与数据大小成正比，那么mae可能会更好。
- mape:平均绝对百分比误差，这是绝对误差与绝对值之间的平均比，这给出了与值相比错误大小的概念。
keras度量标准库包含可以这样调用的mae、mse以及其他误差标准等。
keras.metrics.mean_absolute_error(x_valid,navie_forecast).numpy()
keras.metrics.mean_squared_error(x_valid, naive_forecast).numpy()
移动平均线与差分（Moving average and differencing）
移动平均线（MA）
常见且简单的一种预测方法，即计算移动平均线
移动平均法是时间序列分析中的一种，本质上属于一种高级外推法。通过对历史数据求平均，消除其中的极端值，将预测建立在经过平滑处理后的中间值上。
M
A
1
=
(
P
1
+
P
2
+
⋅
⋅
⋅
+
P
n
)
/
n
MA{_1} = (P{_1}+P{_2}+···+P{_n})/n
MA1=(P1+P2+⋅⋅⋅+Pn)/n
M
A
2
=
(
P
2
+
P
3
+
⋅
⋅
⋅
+
P
n
+
1
)
/
n
MA{_2} = (P{_2}+P{_3}+···+P{_n+1})/n
MA2=(P2+P3+⋅⋅⋅+Pn+1)/n
M
A
3
=
(
P
3
+
P
4
+
⋅
⋅
⋅
+
P
n
+
2
)
/
n
MA{_3} = (P{_3}+P{_4}+···+P{_n+2})/n
MA3=(P3+P4+⋅⋅⋅+Pn+2)/n
…
…
…
其中 M A MA MA代表算术移动平均数.
P n P{_n} Pn代表第n个时间段的数据.
n n n代表需要计算的天数（周期）.
例如1到10十个数字，其平均数便是5.5；而移动则意味着这10个数字的变动.假如第一组是1到10，第二组变动成2到11，第三组又变为3到12，那么，这三组平均数各不相同.而这些不同的平均数的集合，便统称为移动平均数，这些移动平均数构成的线便是移动平均线.
这里黄线是蓝色值平均值的图形，这里一个周期为30天，则它很好地消除了噪音的影响，但是却不能预测出时间序列的季节性与趋势。
差分
当我们拿到数据时，这些数据可能浮动很大，那怎么才能让数据变得平稳呢，其中一种方法就是差分法。
差分法：时间序列在t与t-1时刻的差值。
差分就是将数列中的每一项分别与前一项数做差。
注意得到的差分序列第一个数和原来的第一个数一样（相当于第一个数减0）
差分序列最后比原序列多一个数（相当于0减最后一个数）
预测
构造时间序列
下列代码将为我们构造一个时间序列，包括季节性、趋势以及一些噪声。
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow import keras
def plot_series(time, series, format="-", start=0, end=None):
plt.plot(time[start:end], series[start:end], format)
plt.xlabel("Time")
plt.ylabel("Value")
plt.grid(True)
def trend(time, slope=0):
return slope * time
def seasonal_pattern(season_time):
"""Just an arbitrary pattern, you can change it if you wish"""
return np.where(season_time < 0.4,
np.cos(season_time * 2 * np.pi),
1 / np.exp(3 * season_time))
def seasonality(time, period, amplitude=1, phase=0):
"""Repeats the same pattern at each period"""
season_time = ((time + phase) % period) / period
return amplitude * seasonal_pattern(season_time)
def noise(time, noise_level=1, seed=None):
rnd = np.random.RandomState(seed)
return rnd.randn(len(time)) * noise_level
time = np.arange(4 * 365 + 1, dtype="float32")
baseline = 10
series = trend(time, 0.1)
baseline = 10
amplitude = 40
slope = 0.05
noise_level = 5
# Create the series
series = baseline + trend(time, slope) + seasonality(time, period=365, amplitude=amplitude)
# Update with noise
series += noise(time, noise_level, seed=42)
plt.figure(figsize=(10, 6))
plot_series(time, series)
plt.show()
现在我们获得了时间序列，我们将其分割为训练集和验证集从而进行下一步预测。
split_time = 1000
time_train = time[:split_time]
x_train = series[:split_time]
time_valid = time[split_time:]
x_valid = series[split_time:]
plt.figure(figsize=(10, 6))
plot_series(time_train, x_train)
plt.show()
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plt.show()
时间的前1000作为训练期。
时间在1000之后的部分序列作为验证期，下图当中适当放大，但其季节性与训练期的序列一致。
朴素预测（Naive Forecast）
naive_forecast = series[split_time - 1:-1]
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, naive_forecast)
序列较长，较难观察，现在我们来放大一下验证期的开始部分进行观察。
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid, start=0, end=150)
plot_series(time_valid, naive_forecast, start=1, end=151)
可以观察到在朴素预测中，预测序列落后初始的时间序列一步，这是我们在上面人为设定的。
现在让我们计算验证期内均方误差(mse)和平均绝对误差(mae)：
print(keras.metrics.mean_squared_error(x_valid, naive_forecast).numpy())
print(keras.metrics.mean_absolute_error(x_valid, naive_forecast).numpy())
61.827534
5.937908
这是我们的基准线（baseline），现在让我们尝试使用移动平均线（MA）：
使用移动平均线预测
首先建立移动平均线的函数moving_average_forecast(series, window_size)，它有两个参数，根据输入的时间序列以及窗口数据集大小进行移动平均线的预测，由窗口大小，我们会获得前面数量为窗口大小的序列数据的平均值作为下一个数据的预测，并将预测的数据存储在列表当中。
def moving_average_forecast(series, window_size):
"""Forecasts the mean of the last few values.
If window_size=1, then this is equivalent to naive forecast"""
forecast = []
for time in range(len(series) - window_size):
forecast.append(series[time:time + window_size].mean())
return np.array(forecast)
设置window_size=30，得到预测结果，并绘制在验证期内的预测结果与原实序列数据的图。
moving_avg = moving_average_forecast(series, 30)[split_time - 30:]
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, moving_avg)
计算均方误差(mse)和平均绝对误差(mae)：
print(keras.metrics.mean_squared_error(x_valid, moving_avg).numpy())
print(keras.metrics.mean_absolute_error(x_valid, moving_avg).numpy())
106.674576
7.1424184
我们分析结果，发现这种情况比朴素预测的结果还要糟糕。 考虑到移动平均线无法预测趋势或季节性，因此我们尝试使用差分法去除它们（趋势或季节性）。由于该时间序列的季节性周期为 365 天，我们将从时间 t 的值中减去时间 t – 365 的值进行差分。
diff_series = (series[365:] - series[:-365])
diff_time = time[365:]
plt.figure(figsize=(10, 6))
plot_series(diff_time, diff_series)
plt.show()
得到差分后的图像:
通过观察图中曲线走势，可以看到趋势与季节性在差分过后都已经被消除，接下来我们便采用移动平均线进行预测。
diff_moving_avg = moving_average_forecast(diff_series, 50)[split_time - 365 - 50:]
plt.figure(figsize=(10, 6))
plot_series(time_valid, diff_series[split_time - 365:])
plot_series(time_valid, diff_moving_avg)
plt.show()
由于差分使得原来的时间序列失去了趋势与季节性，所以我们现在通过添加 t – 365 的过去值来恢复它们，作图便得到下列结果。
diff_moving_avg_plus_past = series[split_time - 365:-365] + diff_moving_avg
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, diff_moving_avg_plus_past)
plt.show()
再次计算其均方误差（mse）与平均绝对误差（mae），得到
mse=52.97366
mae=5.839311
误差结果比朴素预测结果要好。然而，在图像当中预测看起来有点太随机了，因为我们只是添加了过去的值，这些值包含噪音很嘈杂。 所以我们对过去的值使用移动平均来消除一些噪音：
diff_moving_avg_plus_smooth_past = moving_average_forecast(series[split_time - 370:-360], 10) + diff_moving_avg
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, diff_moving_avg_plus_smooth_past)
plt.show()
这时的mse=33.452263，mae=4.569442
预测效果变得更优更合理。
三、使用深度学习进行模型的构建
原来，我们只考虑一个合成的季节性的数据集，包括趋势，季节性以及噪声。还学习了分析数据集并且从中作出预测。得到了不错的结果，但我们在其中还并未运用机器学习的内容。
接下来，我们将考虑使用神经网络进行模型的构建与问题的求解。
特征与标签（feature and label）
首先，我们像其他机器学习问题那样需要将数据集分为特征（feature）与标签（label）。在我们进行时间序列预测的情况下，我们的特征是序列当中的很多值，而标签则是这些值的下一个值。我们称前面这些数量的值为我们的特征。
窗口大小（window_size），我们在其中获取数据窗口，训练深度学习模型以预测下一个值。例如，我们获取时间序列数据，例如一次30天，我们用30个值作为特征，下一个值是标签。然后我们训练一个神经网络模型以将30个特征与单个标签进行匹配。
例如我们使用tf.data.Dataset类为我们创造一些数据，我们将设置十个值的范围，当我们输出它们，我们可以看到一系列从0到9的数据。
dataset = tf.data.Dataset.range(10)
for val in dataset:
print(val.numpy())
0
1
2
3
4
5
6
7
8
9
现在，我们使其更加有趣，使用dataset.window通过窗口扩展我们的数据集，它的参数是窗口大小和我们每次想要移动的步数，即偏移量。
这里，我们输出时将窗口大小设置为5且偏移量设置为1。
dataset = tf.data.Dataset.range(10)
dataset = dataset.window(5, shift=1)
for window_dataset in dataset:
for val in window_dataset:
print(val.numpy(), end=" ")
print()
0 1 2 3 4
1 2 3 4 5
2 3 4 5 6
3 4 5 6 7
4 5 6 7 8
5 6 7 8 9
6 7 8 9
7 8 9
8 9
9
观察结果，发现这样窗口末端一旦到了序列数据集的末尾，我们就会得到更少的值，因为后面根本不存在值。所以后面得到6789、789、89等。
因此，我们对窗口进行一些编辑，以便我们可以定期调整数据大小。我们可以在窗口上添加一个名为drop_remainder，如果将其设置为True，便可以通过删除所有剩余部分来截断数据，也就是说它只会给我们五个数据元素的窗口，而其他少于五个数据元素的都不会存在。
dataset = tf.data.Dataset.range(10)
dataset = dataset.window(5, shift=1, drop_remainder=True)
for window_dataset in dataset:
for val in window_dataset:
print(val.numpy(), end=" ")
print()
因此，我们输出它时，它将看起来像下面这样，从01234开始，一直到56789。
0 1 2 3 4
1 2 3 4 5
2 3 4 5 6
3 4 5 6 7
4 5 6 7 8
5 6 7 8 9
现在我们将这些数据放到numpy列表中，这样以后就可以在机器学习中使用这些数据。
dataset = tf.data.Dataset.range(10)
dataset = dataset.window(5, shift=1, drop_remainder=True)
dataset = dataset.flat_map(lambda window: window.batch(5))
for window in dataset:
print(window.numpy())
接下来是将数据分为特征与标签，对于列表中每个元素，所有值都是有意义的，最后一个将作为标签，通过使用map()将列表最后一个值与前面值分开，即lambda window: (window[:-1], window[-1:]）
dataset = tf.data.Dataset.range(10)
dataset = dataset.window(5, shift=1, drop_remainder=True)
dataset = dataset.flat_map(lambda window: window.batch(5))
dataset = dataset.map(lambda window: (window[:-1], window[-1:]))
for x,y in dataset:
print(x.numpy(), y.numpy())
通常，在训练之前要将数据随机打乱洗牌（shuffle），使用 dataset.shuffle()可以实现这一点，我们设置其缓冲区为10，即buffer_sizxe=10，因为这是我们拥有数据的数量。
dataset = tf.data.Dataset.range(10)
dataset = dataset.window(5, shift=1, drop_remainder=True)
dataset = dataset.flat_map(lambda window: window.batch(5))
dataset = dataset.map(lambda window: (window[:-1], window[-1:]))
dataset = dataset.shuffle(buffer_size=10)
for x,y in dataset:
print(x.numpy(), y.numpy())
这时观察到特征与标签已经改组（下述结果不唯一）
[4 5 6 7] [8]
[2 3 4 5] [6]
[5 6 7 8] [9]
[1 2 3 4] [5]
[0 1 2 3] [4]
[3 4 5 6] [7]
最后，我们可以看一下批数据处理，这是通过批处理方法完成的，它需要设置一个size参数，这里设置为2，因此将数据分为两个批进行处理，将它们输出，我们会看到下述结果。
dataset = tf.data.Dataset.range(10)
dataset = dataset.window(5, shift=1, drop_remainder=True)
dataset = dataset.flat_map(lambda window: window.batch(5))
dataset = dataset.map(lambda window: (window[:-1], window[-1:]))
dataset = dataset.shuffle(buffer_size=10)
dataset = dataset.batch(2).prefetch(1)
for x,y in dataset:
print("x = ", x.numpy())
print("y = ", y.numpy())
x = [[3 4 5 6] [1 2 3 4]]
y = [[7] [5]]
x = [[2 3 4 5] [0 1 2 3]]
y = [[6] [4]]
x = [[5 6 7 8] [4 5 6 7]]
y = [[9] [8]]
得到三批数据，每批数据两个数据项，每批数据中相应的数据项都是相对应的。
将窗口数据集输入到神经网络
我们知道了如何创建窗口数据集，接下来，我们将修改代码以供神经网络使用。
我们从此函数开始，该函数将决定如何设置窗口数据集（windowed_dataset），它将接收一个数据系列（series）以及我们想要窗口的大小（window_size），训练时使用的批次大小（batch_size）以及缓冲区大小（shuffle_buffer），它决定如何对将数据进行清洗。
def windowed_dataset(series, window_size, batch_size, shuffle_buffer):
dataset = tf.data.Dataset.from_tensor_slices(series)
dataset = dataset.window(window_size + 1, shift=1, drop_remainder=True)
dataset = dataset.flat_map(lambda window: window.batch(window_size + 1))
dataset = dataset.shuffle(shuffle_buffer).map(lambda window: (window[:-1], window[-1]))
dataset = dataset.batch(batch_size).prefetch(1)
return dataset
第一步是使用 tf.data.Dataset从系列中创建一个数据集，然后使用from_tensor_slices将序列数据传递给它；然后使用数据集的window方法，根据window_size将数据切成适当的窗口，每一个都偏移一个时间，然后drop_remainder设置为True，我们将保持它们每次的大小相同；然后我们将数据展平以使其更易于使用，它将以我们的window_size+1的大小展平为大块；展平后，我们将更简单地将其洗牌，使用洗牌缓冲区将加快处理速度。然后，洗牌后的数据集拆分为xs，它是除最后一个元素外的所有元素，而y是最后一个元素；然后将其分为选定的批次大小并返回。
单层神经网络（Single Neural Network）
现在我们有了一个窗口数据集，我们先看一个超级简单的进行线性回归的神经网络，观察它的准确性，然后从中进行改进。在训练之前，我们需要将数据集拆分为训练期与验证期，这里以1000作为分段时间，前1000作为训练期的数据，1000之后的数据作为验证期的数据。
split_time = 1000
time_train = time[:split_time]
x_train = series[:split_time]
time_valid = time[split_time:]
x_valid = series[split_time:]
我们先从设置的常量进行观察，其中包括我们要训练的窗口大小，数据批量的大小，以及改组之后缓冲区的大小。
window_size = 20
batch_size = 32
shuffle_buffer_size = 1000
然后创建窗口数据集，我们通过采用序列完成这一点，根据提供的各个常量，提供格式化的数据集。
然后建立一个神经网络的dense层，输入大小为window_size，并将该层只有一个神经元的神经网络给予一个变量名为l0，因为之后要输出它的权重，这样会简单的多。
然后建立单个变量的顺序模型。
dataset = windowed_dataset(series, window_size, batch_size, shuffle_buffer_size)
l0 = tf.keras.layers.Dense(1, input_shape=[window_size])
model = tf.keras.models.Sequential([l0])
通过编译，设置均方误差mse作为损失函数，以随机梯度下降作为优化方法（这种优化方法有着可以自己设定参数等优点），最后进行拟合，使用以及格式化的数据集dataset，设置拟合的次数100次。
model.compile(loss="mse", optimizer=tf.keras.optimizers.SGD(lr=1e-6, momentum=0.9))
model.fit(dataset,epochs=100,verbose=0)
完成训练后，检查神经网络的不同权重，对变量l0命名就是为了实现下述代码
print("Layer weights {}".format(l0.get_weights()))
Layer weights [array([[-0.05529308],
[ 0.00947868],
[ 0.05041581],
[-0.00872855],
[-0.02697379],
[ 0.04566243],
[-0.00740558],
[ 0.00707938],
[-0.00853621],
[ 0.01960586],
[-0.03908759],
[-0.0327211 ],
[ 0.03778784],
[ 0.07210432],
[-0.05659156],
[ 0.09618359],
[-0.01626248],
[ 0.22845924],
[ 0.17617637],
[ 0.5097049 ]], dtype=float32), array([0.01695044], dtype=float32)]
可以观察到，第一个数组中有20个值，而第二个数组中只有一个值，因为神经网络在训练之后已经了解了线性回归，所以20个值已经可以作为20个值的权重
w
n
w{_n}
wn，另一个则为偏差值
b
b
b。
Y
=
w
0
x
0
+
w
1
x
1
+
w
2
x
2
+
⋅
⋅
⋅
+
w
1
9
x
1
9
+
b
Y=w{_0}x{_0}+w{_1}x{_1}+w{_2}x{_2}+···+w{_1}{_9}x{_1}{_9}+b
Y=w0x0+w1x1+w2x2+⋅⋅⋅+w19x19+b
我们考虑一下输入窗口宽度的值，为
x
0
x{_0}
x0、
x
1
x{_1}
x1一直到
x
1
9
x{_1}{_9}
x19，但我们需要清楚一点，那并不能作为水平轴上的值，而是时间序列的值在水平轴上那个点，因此在时间
t
0
t{_0}
t0的那个值
现在，让我们看一下这些值是在将当前值称为 x 0 x{_0} x0的前20步，然后依次类推。
然后我们的预测就是权重乘上x的值然后再进行加上偏差，例如我拿序列的20项并且输出它们，我们可以看到20个值，如果我想预测它们，我们将这些值传递到模型当中，然后使用np.newaxis改变其输入尺寸以供模型使用（np.newaxis将会额外生成一个维度）
print(series[1:21])
model.predict(series[1:21][np.newaxis])
[49.35275 53.314735 57.711823 48.934444 48.931244 57.982895 53.897125
47.67393 52.68371 47.591717 47.506374 50.959415 40.086178 40.919415
46.612473 44.228207 50.720642 44.454983 41.76799 55.980938]
array([[50.9769]], dtype=float32)
上述结果的上面的数据是我们选用的20个值，下面的数据是我们返回的预测值，也就是说，我们预测的模型看到像上面这样20个值会预测得到50.9769这样一个值。因此，如果我们要绘制对每个点的预测，我们可以构造一个放置预测数据的空列表，然后开始遍历序列，以时间与窗口大小作为切片，np.newaxis将序列尺寸改为可供模型预测的尺寸，每次遍历预测之后，将预测值添加到列表中。
我们在验证期进行验证，选取split_time-window_size预测到的值与实际相对照，得到下图。
forecast = []
for time in range(len(series) - window_size):
forecast.append(model.predict(series[time:time + window_size][np.newaxis]))
forecast = forecast[split_time-window_size:]
results = np.array(forecast)[:, 0, 0]
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, results)
橙色线为预测线，蓝色线为原始序列数据。
计算其平均绝对误差（mae），从而观察它的预测效果
tf.keras.metrics.mean_absolute_error(x_valid, results).numpy()
使用深层神经网络（Using Deep Neural Network）
上面我们只是考虑了使用单个神经元进行线性回归的简单方法，而接下来我们考虑使用深层的神经网络，看是否能得到一个寻来你效果更佳的模型。
首先构造神经网络，这里一共构造了三层，第一层为Dense层，设置其包含10个神经元，输入尺寸为窗口数据大小，第二层相同，第三层包含单个神经元，它得到的值便是我们通过预测所要获得的值，这三层神经网络均使用relu函数作为机会哦函数进行激活。
然后编译，采用均方误差作为损失函数，采用随机梯度下降法作为优化方法，最后进行拟合训练100次，得到训练好的模型。
其他内容与上面内容相同，输出预测结果，得到平均绝对误差（mae）为5.2200565
dataset = windowed_dataset(x_train, window_size, batch_size, shuffle_buffer_size)
model = tf.keras.models.Sequential([
tf.keras.layers.Dense(10, input_shape=[window_size], activation="relu"),
tf.keras.layers.Dense(10, activation="relu"),
tf.keras.layers.Dense(1)
])
model.compile(loss="mse", optimizer=tf.keras.optimizers.SGD(lr=1e-6, momentum=0.9))
model.fit(dataset,epochs=100,verbose=0)
forecast = []
for time in range(len(series) - window_size):
forecast.append(model.predict(series[time:time + window_size][np.newaxis]))
forecast = forecast[split_time-window_size:]
results = np.array(forecast)[:, 0, 0]
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, results)
tf.keras.metrics.mean_absolute_error(x_valid, results).numpy()
5.2200565
为了使我们的模型训练效果更好，下面我们使用了callbacks，从而为我们找到一个最佳的学习率（learning_rate）
dataset = windowed_dataset(x_train, window_size, batch_size, shuffle_buffer_size)
model = tf.keras.models.Sequential([
tf.keras.layers.Dense(10, input_shape=[window_size], activation="relu"),
tf.keras.layers.Dense(10, activation="relu"),
tf.keras.layers.Dense(1)
])
lr_schedule = tf.keras.callbacks.LearningRateScheduler(
lambda epoch: 1e-8 * 10**(epoch / 20))
optimizer = tf.keras.optimizers.SGD(lr=1e-8, momentum=0.9)
model.compile(loss="mse", optimizer=optimizer)
history = model.fit(dataset, epochs=100, callbacks=[lr_schedule], verbose=0)
lrs = 1e-8 * (10 ** (np.arange(100) / 20))
plt.semilogx(lrs, history.history["loss"])
plt.axis([1e-8, 1e-3, 0, 300])
[1e-08, 0.001, 0, 300]
通过观察学习率与损失函数值的关系图，x轴显示学习率，y轴显示的是损失函数值，我们可以看到学习率在8e-6左右时，损失函数值最小。
window_size = 30
dataset = windowed_dataset(x_train, window_size, batch_size, shuffle_buffer_size)
model = tf.keras.models.Sequential([
tf.keras.layers.Dense(10, activation="relu", input_shape=[window_size]),
tf.keras.layers.Dense(10, activation="relu"),
tf.keras.layers.Dense(1)
])
optimizer = tf.keras.optimizers.SGD(lr=8e-6, momentum=0.9)
model.compile(loss="mse", optimizer=optimizer)
history = model.fit(dataset, epochs=500, verbose=0)
loss = history.history['loss']
epochs = range(len(loss))
plt.plot(epochs, loss, 'b', label='Training Loss')
plt.show()
我们就可以看到损失函数值随时间的变化情况，可以看到其急剧下降并且变平，可能只是训练了10次左右，损失值就变得很平不能观察。
所以我们通过删掉前十次的训练的损失值进行绘制图表。
# Plot all but the first 10
loss = history.history['loss']
epochs = range(10, len(loss))
plot_loss = loss[10:]
print(plot_loss)
plt.plot(epochs, plot_loss, 'b', label='Training Loss')
plt.show()
我们可以看到，经历了500次训练之后，损失值仍然在减小，说明神经网络确实学习的很好。
将预测结果与原始数据对比可得下图。
forecast = []
for time in range(len(series) - window_size):
forecast.append(model.predict(series[time:time + window_size][np.newaxis]))
forecast = forecast[split_time-window_size:]
results = np.array(forecast)[:, 0, 0]
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, results)
计算得到其mae也变得更小，说明模型得到了更好的效果。
tf.keras.metrics.mean_absolute_error(x_valid, results).numpy()
4.8788633
RNN和LSTM的使用
原来，我们已经了解了如何使用DNN进行神经网络的构建，从一些简单的分析技术开始，转化为使用机器学习做一个简单的回归，并且使用DNN做了一些简单的调整，得到了效果更好的模型，下面，我们将使用RNN进行预测。
循环神经网络（Recurrent Neural Network）与其他神经网络不同的是它包括了循环层，这是为了能够顺利地处理输入的序列，RNN非常灵活，我们可以用它预测文本等，这里我们只用它处理时间序列。我们下面要构建一个包括两个循环层以及一个致密层的RNN神经网络。使用RNN，可以进行分批填充，它也将输出一批预测，想DNN得到的结果那样。有一个不同之处就是RNN完整的输入形状是三维的，shape=[batch_size,#time_steps,#dims]，第一个维度是批量大小，第二个维度是时间步，第三个是维度。
例如，单变量时间序列维度值将为1，对于多变量，它将更多。
接下来，让我们研究一下RNN的工作原理看看它们到底如何工作。下图看起来有很多细胞，但实际上只有一个并反复用于计算输出，在此图中看起来有很多但使用的一直只有一个，在每个时间步，存储单元将获取该步骤的输入值，例如在0时间它为0，输入0状态并且输入X0得到Y0，并传递向量H0用于下一步计算，H0与X1一起馈入单元格输出Y1，以此类推。
RNN输入形状
输入是三维的，例如我们有30个时间步，我们按四个大小批量处理它们，它的维度为1，这样我们的输入形状便为4301.存储单元得到的将是下图这样一个4*1矩阵，单元格也将接受上一步的状态矩阵，在第一步中，它将为0，对于后续的是存储单元的输出，但除了状态矩阵，我们还会得到一个Y值。
我们可以知道，它的输出是4*3矩阵，因为它的输入是四个数据，而每个单元格中包含三个神经元，因此该层的完整输出是三维的。状态矩阵H只是输出矩阵Y的一个副本，H0是Y0的一个副本，而H1是Y1的一个副本，依次类推。
而有时候我们想要输入一个序列，不需要其持续输出，而是得到我们最终想要获得的向量，这通常称为向量RNN的序列，但实际上要做的就是忽略除了最后一个输出的所有输出。在Tensorflow以及keras当中这是默认的，因此我们想要输出序列时还应该进行设置。
输出一个序列
model = tf.keras.models.Sequential([
tf.keras.layers.SimpleRNN(40, return_sequences=True,input_shape=[None,1]),
tf.keras.layers.SimpleRNN(40),
tf.keras.layers.Dense(1),
])
第一层设置了return_sequences=True，则我们将输出一系列数据，下一层没有设置return_sequences=True，则直接输出最后一个dense层。
看一下input_shape，第一个维度默认为批量大小，因为可以有任意大小，所以不用设置它，第二个维度为时间步数，我们可以将其设置为None，这意味着我们可以处理任意长度的时间序列，最后一个维度为1，因为我们使用的是时间序列的单位变化。
若第二层也设置return_sequences=True，最后的结果将返回一个长度相同的序列。
Lambda Layer
使用Lambda Layer能够帮助我们拓展Tensorflow的keras功能，我们可以在模型定义本身中做到这一点，第一个Lambda Layer将为我们提高尺寸，将数据扩展一维。input_shape=[None]意味这可以处理任意大小的序列。
同样的，如果我们将输出放大100倍，有助于我们后续的训练过程。
RNN默认的激活函数为双曲正切函数tanh，这将输出-1到1之间的值。
model = tf.keras.models.Sequential([
tf.keras.layers.Lambda(lambda x: tf.expand_dims(x, axis=-1),
input_shape=[None]),
tf.keras.layers.SimpleRNN(40, return_sequences=True),
tf.keras.layers.SimpleRNN(40),
tf.keras.layers.Dense(1),
tf.keras.layers.Lambda(lambda x: x * 100.0)
])
调整学习率
使用callbacks，训练100次，绘制出图形。
train_set = windowed_dataset(x_train, window_size, batch_size=128, shuffle_buffer=shuffle_buffer_size)
model = tf.keras.models.Sequential([
tf.keras.layers.Lambda(lambda x: tf.expand_dims(x, axis=-1),
input_shape=[None]),
tf.keras.layers.SimpleRNN(40, return_sequences=True),
tf.keras.layers.SimpleRNN(40),
tf.keras.layers.Dense(1),
tf.keras.layers.Lambda(lambda x: x * 100.0)
])
lr_schedule = tf.keras.callbacks.LearningRateScheduler(
lambda epoch: 1e-8 * 10**(epoch / 20))
optimizer = tf.keras.optimizers.SGD(lr=1e-8, momentum=0.9)
model.compile(loss=tf.keras.losses.Huber(),
optimizer=optimizer,
metrics=["mae"])
history = model.fit(train_set, epochs=100, callbacks=[lr_schedule])
plt.semilogx(history.history["lr"], history.history["loss"])
plt.axis([1e-8, 1e-4, 0, 30])
通过观察损失函数值与不同学习率对应关系，可以看到造学习率等于5e-6左右损失值最小。
RNN
对模型使用调整好的学习率训练400次。
tf.keras.backend.clear_session()
tf.random.set_seed(51)
np.random.seed(51)
dataset = windowed_dataset(x_train, window_size, batch_size=128, shuffle_buffer=shuffle_buffer_size)
model = tf.keras.models.Sequential([
tf.keras.layers.Lambda(lambda x: tf.expand_dims(x, axis=-1),
input_shape=[None]),
tf.keras.layers.SimpleRNN(40, return_sequences=True),
tf.keras.layers.SimpleRNN(40),
tf.keras.layers.Dense(1),
tf.keras.layers.Lambda(lambda x: x * 100.0)
])
optimizer = tf.keras.optimizers.SGD(lr=5e-5, momentum=0.9)
model.compile(loss=tf.keras.losses.Huber(),
optimizer=optimizer,
metrics=["mae"])
history = model.fit(dataset,epochs=400)
训练结束后我们便可以去验证验证集并且将其结果进行绘制。
forecast=[]
for time in range(len(series) - window_size):
forecast.append(model.predict(series[time:time + window_size][np.newaxis]))
forecast = forecast[split_time-window_size:]
results = np.array(forecast)[:, 0, 0]
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, results)
观察预测结果与原始数据可以看出，除了波峰位置预测不太尽如人意，其他部分的预测也都还不错。
我们输出其mae发现值为6.9596496，结果还可以。
tf.keras.metrics.mean_absolute_error(x_valid, results).numpy()
6.9596496
我们再来看一下损失函数值loss与mae和训练次数的关系图。
loss与mae在50次之后仍在下降，但在400次左右就变得不再非常稳定，所以我们在400次将其终止。
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
#-----------------------------------------------------------
# Retrieve a list of list results on training and test data
# sets for each training epoch
#-----------------------------------------------------------
mae=history.history['mae']
loss=history.history['loss']
epochs=range(len(loss)) # Get number of epochs
#------------------------------------------------
# Plot MAE and Loss
#------------------------------------------------
plt.plot(epochs, mae, 'r')
plt.plot(epochs, loss, 'b')
plt.title('MAE and Loss')
plt.xlabel("Epochs")
plt.ylabel("Accuracy")
plt.legend(["MAE", "Loss"])
plt.figure()
epochs_zoom = epochs[200:]
mae_zoom = mae[200:]
loss_zoom = loss[200:]
#------------------------------------------------
# Plot Zoomed MAE and Loss
#------------------------------------------------
plt.plot(epochs_zoom, mae_zoom, 'r')
plt.plot(epochs_zoom, loss_zoom, 'b')
plt.title('MAE and Loss')
plt.xlabel("Epochs")
plt.ylabel("Accuracy")
plt.legend(["MAE", "Loss"])
plt.figure()
LSTM
神经网络会受到短时记忆的影响。如果一条序列足够长，那它们将很难将信息从较早的时间步传送到后面的时间步。 因此，如果你正在尝试处理一段长序列进行预测，RNN 可能从一开始就会遗漏重要信息。
在之前，我们使用RNN进行了一些预测，结果不错，但是我们还是需要进行一些改善，因为我们在预测当中出现了比较奇怪的波峰，尝试使用不同的超参数，我们将看到一些改进，但更好的方法是使用LSTM去代替RNN来查看影响。
在RNN当中，它的特点会成为影响后续计算的一个因素，它的影响会随着时间推移大大减小，LSTM是将在神经网络当中一直保持状态，状态从一个单元传递到另一个单元，可以被更好地进行维护，这说明早先的数据对总体预测的影响比RNN要高。状态也可以是向前或者是向后的，对文本来说这非常重要，因为文本需要结合上下文进行分析，对数字序列的影响可能较小，关于LSTM内容较多，建议从其他资料进行了解其原理，这里仅作简单介绍。
下面分析代码，并且查看LSTM及其工作方式.
首先是tf.keras.backend.clear_session()，这将会清除所有的内部变量，这样就可以放心地进行实验而不用担心影响到其他的部分。
在使用Lambda层修改完尺寸之后，我们添加了具有32个单元格的LSTM层，我们也对其使用双向来观察其对预测产生的影响。
输出神经元为我们提供预测值。
tf.keras.backend.clear_session()
tf.random.set_seed(51)
np.random.seed(51)
tf.keras.backend.clear_session()
dataset = windowed_dataset(x_train, window_size, batch_size, shuffle_buffer_size)
model = tf.keras.models.Sequential([
tf.keras.layers.Lambda(lambda x: tf.expand_dims(x, axis=-1),
input_shape=[None]),
tf.keras.layers.Bidirectional(tf.keras.layers.LSTM(32, return_sequences=True)),
tf.keras.layers.Bidirectional(tf.keras.layers.LSTM(32)),
tf.keras.layers.Dense(1),
tf.keras.layers.Lambda(lambda x: x * 100.0)
])
lr_schedule = tf.keras.callbacks.LearningRateScheduler(
lambda epoch: 1e-8 * 10**(epoch / 20))
optimizer = tf.keras.optimizers.SGD(lr=1e-8, momentum=0.9)
model.compile(loss=tf.keras.losses.Huber(),
optimizer=optimizer,
metrics=["mae"])
history = model.fit(dataset, epochs=100, callbacks=[lr_schedule])
forecast = []
results = []
for time in range(len(series) - window_size):
forecast.append(model.predict(series[time:time + window_size][np.newaxis]))
forecast = forecast[split_time-window_size:]
results = np.array(forecast)[:, 0, 0]
plt.figure(figsize=(10, 6))
plot_series(time_valid, x_valid)
plot_series(time_valid, results)
tf.keras.metrics.mean_absolute_error(x_valid, results).numpy()
更新中…
6622
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
