Machine-Learning-Mastery-深度学习时间序列教程-十-

Machine Learning Mastery 深度学习时间序列教程（十）

原文：Machine Learning Mastery

协议：CC BY-NC-SA 4.0

如何在 Python 中为时间序列预测网格搜索三次指数平滑

原文： machinelearningmastery.com/how-to-grid-search-triple-exponential-smoothing-for-time-series-forecasting-in-python/

指数平滑是单变量数据的时间序列预测方法，可以扩展为支持具有系统趋势或季节性成分的数据。

通常的做法是使用优化过程来查找模型超参数，这些参数导致指数平滑模型具有给定时间序列数据集的最佳表现。此实践仅适用于模型用于描述水平，趋势和季节性的指数结构的系数。

还可以自动优化指数平滑模型的其他超参数，例如是否对趋势和季节性分量建模，如果是，是否使用加法或乘法方法对它们进行建模。

在本教程中，您将了解如何开发一个框架，用于网格搜索所有指数平滑模型超参数，以进行单变量时间序列预测。

完成本教程后，您将了解：

• 如何使用前向验证从零开始开发网格搜索 ETS 模型的框架。
• 如何为女性出生日常时间序列数据网格搜索 ETS 模型超参数。
• 如何针对洗发水销售，汽车销售和温度的月度时间序列数据网格搜索 ETS 模型超参数。

让我们开始吧。

• Oct8 / 2018 ：更新了 ETS 模型的拟合，以使用 NumPy 数组修复乘法趋势/季节性问题（感谢 Amit Amola）。

如何在 Python 中进行时间序列预测的网格搜索三次指数平滑 照片由 john mcsporran 拍摄，保留一些权利。

教程概述

本教程分为六个部分;他们是：

1. 时间序列预测的指数平滑
2. 开发网格搜索框架
3. 案例研究 1：没有趋势或季节性
4. 案例研究 2：趋势
5. 案例研究 3：季节性
6. 案例研究 4：趋势和季节性

时间序列预测的指数平滑

指数平滑是单变量数据的时间序列预测方法。

像 Box-Jenkins ARIMA 系列方法这样的时间序列方法开发了一种模型，其中预测是近期过去观察或滞后的加权线性和。

指数平滑预测方法的类似之处在于预测是过去观察的加权和，但模型明确地使用指数减小的权重用于过去的观察。

具体而言，过去的观察以几何减小的比率加权。

使用指数平滑方法产生的预测是过去观测的加权平均值，随着观测结果的变化，权重呈指数衰减。换句话说，观察越近，相关重量越高。

• 第 171 页，预测：原则和实践，2013。

指数平滑方法可以被视为对等，并且是流行的 Box-Jenkins ARIMA 类时间序列预测方法的替代方法。

总的来说，这些方法有时被称为 ETS 模型，参考 _ 错误 ， 趋势 _ 和 _ 季节性 _ 的显式建模。

指数平滑有三种类型;他们是：

• 单指数平滑或 SES，用于没有趋势或季节性的单变量数据。
• 双指数平滑用于支持趋势的单变量数据。
• 三重指数平滑，或 Holt-Winters 指数平滑，支持趋势和季节性。

三指数平滑模型通过趋势性质（加法，乘法或无）的性质和季节性的性质（加法，乘法或无）来表示单指数和双指数平滑，以及任何阻尼趋势。

开发网格搜索框架

在本节中，我们将为给定的单变量时间序列预测问题开发一个网格搜索指数平滑模型超参数的框架。

我们将使用 statsmodels 库提供的 Holt-Winters 指数平滑的实现。

该模型具有超参数，可控制为系列，趋势和季节性执行的指数的性质，具体为：

• smoothing_level （alpha）：该级别的平滑系数。
• smoothing_slope （beta）：趋势的平滑系数。
• smoothing_seasonal （gamma）：季节性成分的平滑系数。
• damping_slope （phi）：阻尼趋势的系数。

在定义模型时，可以指定所有这四个超参数。如果未指定它们，库将自动调整模型并找到这些超参数的最佳值（例如 optimized = True ）。

还有其他超参数，模型不会自动调整您可能想要指定的;他们是：

• 趋势：趋势分量的类型，作为加法的“_ 加 _”或乘法的“mul”。可以通过将趋势设置为“无”来禁用对趋势建模。
• 阻尼：趋势分量是否应该被阻尼，无论是真还是假。
• 季节性：季节性成分的类型，为“_ 添加 _”为添加剂或“mul”为乘法。可以通过将季节性组件设置为“无”来禁用它。
• seasonal_periods ：季节性时间段内的时间步数，例如在一年一度的季节性结构中 12 个月 12 个月。
• use_boxcox ：是否执行系列的幂变换（True / False）或指定变换的 lambda。

如果您对问题了解得足以指定其中一个或多个参数，则应指定它们。如果没有，您可以尝试网格搜索这些参数。

我们可以通过定义一个适合具有给定配置的模型的函数来开始，并进行一步预测。

下面的exp_smoothing_forecast()实现了这种行为。

该函数采用连续先前观察的数组或列表以及用于配置模型的配置参数列表。

配置参数依次为：趋势类型，阻尼类型，季节性类型，季节周期，是否使用 Box-Cox 变换，以及在拟合模型时是否消除偏差。

# one-step Holt Winter's Exponential Smoothing forecast
def exp_smoothing_forecast(history, config):
	t,d,s,p,b,r = config
	# define model model
	history = array(history)
	model = ExponentialSmoothing(history, trend=t, damped=d, seasonal=s, seasonal_periods=p)
	# fit model
	model_fit = model.fit(optimized=True, use_boxcox=b, remove_bias=r)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

接下来，我们需要建立一些函数，通过前向验证重复拟合和评估模型，包括将数据集拆分为训练和测试集并评估一步预测。

我们可以使用给定指定大小的分割的切片来分割列表或 NumPy 数据数组，例如，从测试集中的数据中使用的时间步数。

下面的train_test_split()函数为提供的数据集和要在测试集中使用的指定数量的时间步骤实现此功能。

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

在对测试数据集中的每个步骤做出预测之后，需要将它们与测试集进行比较以计算错误分数。

时间序列预测有许多流行的错误分数。在这种情况下，我们将使用均方根误差（RMSE），但您可以将其更改为您的首选度量，例如 MAPE，MAE 等

下面的measure_rmse()函数将根据实际（测试集）和预测值列表计算 RMSE。

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

我们现在可以实现前向验证方案。这是评估尊重观测时间顺序的时间序列预测模型的标准方法。

首先，使用train_test_split()函数将提供的单变量时间序列数据集分成训练集和测试集。然后枚举测试集中的观察数。对于每一个，我们在所有历史记录中拟合模型并进行一步预测。然后将对时间步骤的真实观察添加到历史中，并重复该过程。调用exp_smoothing_forecast()函数以适合模型并做出预测。最后，通过调用measure_rmse()函数，将所有一步预测与实际测试集进行比较，计算错误分数。

下面的walk_forward_validation()函数实现了这一点，采用了单变量时间序列，在测试集中使用的一些时间步骤，以及一组模型配置。

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = exp_smoothing_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

如果您对进行多步预测感兴趣，可以在exp_smoothing_forecast()函数中更改predict()的调用，并更改 _ 中的错误计算 measure_rmse（）_ 功能。

我们可以使用不同的模型配置列表重复调用 walk_forward_validation（）。

一个可能的问题是，可能不会为模型调用模型配置的某些组合，并且会抛出异常，例如，指定数据中季节性结构的一些但不是所有方面。

此外，某些型号还可能会对某些数据发出警告，例如：来自 statsmodels 库调用的线性代数库。

我们可以在网格搜索期间捕获异常并忽略警告，方法是将所有调用包含在walk_forward_validation()中，并使用 try-except 和 block 来忽略警告。我们还可以添加调试支持来禁用这些保护，以防我们想要查看实际情况。最后，如果确实发生错误，我们可以返回 _ 无 _ 结果;否则，我们可以打印一些关于评估的每个模型的技能的信息。当评估大量模型时，这很有用。

下面的score_model()函数实现了这个并返回（键和结果）的元组，其中键是测试模型配置的字符串版本。

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

接下来，我们需要一个循环来测试不同模型配置的列表。

这是驱动网格搜索过程的主要功能，并将为每个模型配置调用score_model()函数。

通过并行评估模型配置，我们可以大大加快网格搜索过程。一种方法是使用 Joblib 库。

我们可以定义一个Parallel对象，其中包含要使用的核心数，并将其设置为硬件中检测到的 CPU 核心数。

executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')

然后我们可以创建一个并行执行的任务列表，这将是对我们拥有的每个模型配置的score_model()函数的一次调用。

tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)

最后，我们可以使用Parallel对象并行执行任务列表。

scores = executor(tasks)

而已。

我们还可以提供评估所有模型配置的非并行版本，以防我们想要调试某些内容。

scores = [score_model(data, n_test, cfg) for cfg in cfg_list]

评估配置列表的结果将是元组列表，每个元组都有一个名称，该名称总结了特定的模型配置，并且使用该配置评估的模型的错误为 RMSE，如果出现错误则为 None。

我们可以使用“无”过滤掉所有分数。

scores = [r for r in scores if r[1] != None]

然后我们可以按照升序排列列表中的所有元组（最好是第一个），然后返回此分数列表以供审阅。

给定单变量时间序列数据集，模型配置列表（列表列表）以及在测试集中使用的时间步数，下面的grid_search()函数实现此行为。可选的并行参数允许对所有内核的模型进行开启或关闭调整，默认情况下处于打开状态。

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

我们差不多完成了。

剩下要做的唯一事情是定义模型配置列表以尝试数据集。

我们可以一般地定义它。我们可能想要指定的唯一参数是系列中季节性组件的周期性（如果存在）。默认情况下，我们假设没有季节性组件。

下面的exp_smoothing_configs()函数将创建要评估的模型配置列表。

可以指定季节性时段的可选列表，您甚至可以更改该功能以指定您可能了解的有关时间序列的其他元素。

从理论上讲，有 72 种可能的模型配置需要评估，但在实践中，许多模型配置无效并会导致我们将陷入和忽略的错误。

# create a set of exponential smoothing configs to try
def exp_smoothing_configs(seasonal=[None]):
	models = list()
	# define config lists
	t_params = ['add', 'mul', None]
	d_params = [True, False]
	s_params = ['add', 'mul', None]
	p_params = seasonal
	b_params = [True, False]
	r_params = [True, False]
	# create config instances
	for t in t_params:
		for d in d_params:
			for s in s_params:
				for p in p_params:
					for b in b_params:
						for r in r_params:
							cfg = [t,d,s,p,b,r]
							models.append(cfg)
	return models

我们现在有一个网格搜索三重指数平滑模型超参数的框架，通过一步前进验证。

它是通用的，适用于作为列表或 NumPy 数组提供的任何内存中单变量时间序列。

我们可以通过在人为设计的 10 步数据集上进行测试来确保所有部分协同工作。

下面列出了完整的示例。

# grid search holt winter's exponential smoothing
from math import sqrt
from multiprocessing import cpu_count
from joblib import Parallel
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.metrics import mean_squared_error
from numpy import array

# one-step Holt Winter’s Exponential Smoothing forecast
def exp_smoothing_forecast(history, config):
	t,d,s,p,b,r = config
	# define model
	history = array(history)
	model = ExponentialSmoothing(history, trend=t, damped=d, seasonal=s, seasonal_periods=p)
	# fit model
	model_fit = model.fit(optimized=True, use_boxcox=b, remove_bias=r)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = exp_smoothing_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

# create a set of exponential smoothing configs to try
def exp_smoothing_configs(seasonal=[None]):
	models = list()
	# define config lists
	t_params = ['add', 'mul', None]
	d_params = [True, False]
	s_params = ['add', 'mul', None]
	p_params = seasonal
	b_params = [True, False]
	r_params = [True, False]
	# create config instances
	for t in t_params:
		for d in d_params:
			for s in s_params:
				for p in p_params:
					for b in b_params:
						for r in r_params:
							cfg = [t,d,s,p,b,r]
							models.append(cfg)
	return models

if`_name_`== '__main__':
	# define dataset
	data = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0]
	print(data)
	# data split
	n_test = 4
	# model configs
	cfg_list = exp_smoothing_configs()
	# grid search
	scores = grid_search(data, cfg_list, n_test)
	print('done')
	# list top 3 configs
	for cfg, error in scores[:3]:
		print(cfg, error)

首先运行该示例打印设计的时间序列数据集。

接下来，在评估模型配置及其错误时报告它们。

最后，报告前三种配置的配置和错误。

[10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0]

 > Model[[None, False, None, None, True, True]] 1.380
 > Model[[None, False, None, None, True, False]] 10.000
 > Model[[None, False, None, None, False, True]] 2.563
 > Model[[None, False, None, None, False, False]] 10.000
done

[None, False, None, None, True, True] 1.379824445857423
[None, False, None, None, False, True] 2.5628662672606612
[None, False, None, None, False, False] 10.0

我们不报告模型本身优化的模型参数。假设您可以通过指定更广泛的超参数来再次获得相同的结果，并允许库找到相同的内部参数。

您可以通过重新配置具有相同配置的独立模型并在模型拟合上打印'params'属性的内容来访问这些内部参数;例如：

print(model_fit.params)

现在我们有了一个强大的网格搜索框架来搜索 ETS 模型超参数，让我们在一套标准的单变量时间序列数据集上进行测试。

选择数据集用于演示目的;我并不是说 ETS 模型是每个数据集的最佳方法，在某些情况下，SARIMA 或其他东西可能更合适。

案例研究 1：没有趋势或季节性

“每日女性分娩”数据集总结了 1959 年美国加利福尼亚州每日女性总分娩数。

数据集没有明显的趋势或季节性成分。

每日女性出生数据集的线图

您可以从 DataMarket 了解有关数据集的更多信息。

直接从这里下载数据集：

• 每日总数 - 女性分娩.sv

在当前工作目录中使用文件名“ daily-total-female-births.csv ”保存文件。

我们可以使用函数read_csv()将此数据集作为 Pandas 系列加载。

series = read_csv('daily-total-female-births.csv', header=0, index_col=0)

数据集有一年或 365 个观测值。我们将使用前 200 个进行训练，将剩余的 165 个作为测试集。

下面列出了搜索每日女性单变量时间序列预测问题的完整示例网格。

# grid search ets models for daily female births
from math import sqrt
from multiprocessing import cpu_count
from joblib import Parallel
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.metrics import mean_squared_error
from pandas import read_csv
from numpy import array

# one-step Holt Winter’s Exponential Smoothing forecast
def exp_smoothing_forecast(history, config):
	t,d,s,p,b,r = config
	# define model
	history = array(history)
	model = ExponentialSmoothing(history, trend=t, damped=d, seasonal=s, seasonal_periods=p)
	# fit model
	model_fit = model.fit(optimized=True, use_boxcox=b, remove_bias=r)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = exp_smoothing_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

# create a set of exponential smoothing configs to try
def exp_smoothing_configs(seasonal=[None]):
	models = list()
	# define config lists
	t_params = ['add', 'mul', None]
	d_params = [True, False]
	s_params = ['add', 'mul', None]
	p_params = seasonal
	b_params = [True, False]
	r_params = [True, False]
	# create config instances
	for t in t_params:
		for d in d_params:
			for s in s_params:
				for p in p_params:
					for b in b_params:
						for r in r_params:
							cfg = [t,d,s,p,b,r]
							models.append(cfg)
	return models

if`_name_`== '__main__':
	# load dataset
	series = read_csv('daily-total-female-births.csv', header=0, index_col=0)
	data = series.values
	# data split
	n_test = 165
	# model configs
	cfg_list = exp_smoothing_configs()
	# grid search
	scores = grid_search(data[:,0], cfg_list, n_test)
	print('done')
	# list top 3 configs
	for cfg, error in scores[:3]:
		print(cfg, error)

运行该示例可能需要几分钟，因为在现代硬件上安装每个 ETS 模型大约需要一分钟。

在评估模型时打印模型配置和 RMSE 在运行结束时报告前三个模型配置及其错误。

我们可以看到最好的结果是大约 6.96 个出生的 RMSE，具有以下配置：

• 趋势：乘法
• 阻尼：错误
• 季节性：无
• 季节性时期：无
• Box-Cox 变换：是的
• 删除偏差：是的

令人惊讶的是，假设乘法趋势的模型比不具有乘法趋势的模型表现得更好。

除非我们抛弃假设和网格搜索模型，否则我们不会知道情况就是这样。

 > Model[['add', False, None, None, True, True]] 7.081
 > Model[['add', False, None, None, True, False]] 7.113
 > Model[['add', False, None, None, False, True]] 7.112
 > Model[['add', False, None, None, False, False]] 7.115
 > Model[['add', True, None, None, True, True]] 7.118
 > Model[['add', True, None, None, True, False]] 7.170
 > Model[['add', True, None, None, False, True]] 7.113
 > Model[['add', True, None, None, False, False]] 7.126
 > Model[['mul', True, None, None, True, True]] 7.118
 > Model[['mul', True, None, None, True, False]] 7.170
 > Model[['mul', True, None, None, False, True]] 7.113
 > Model[['mul', True, None, None, False, False]] 7.126
 > Model[['mul', False, None, None, True, True]] 6.961
 > Model[['mul', False, None, None, True, False]] 6.985
 > Model[[None, False, None, None, True, True]] 7.169
 > Model[[None, False, None, None, True, False]] 7.212
 > Model[[None, False, None, None, False, True]] 7.117
 > Model[[None, False, None, None, False, False]] 7.126
done

['mul', False, None, None, True, True] 6.960703917145126
['mul', False, None, None, True, False] 6.984513598720297
['add', False, None, None, True, True] 7.081359856193836

案例研究 2：趋势

“洗发水”数据集总结了三年内洗发水的月销售额。

数据集包含明显的趋势，但没有明显的季节性成分。

月度洗发水销售数据集的线图

您可以从 DataMarket 了解有关数据集的更多信息。

直接从这里下载数据集：

• shampoo.csv

在当前工作目录中使用文件名“shampoo.csv”保存文件。

我们可以使用函数read_csv()将此数据集作为 Pandas 系列加载。

# parse dates
def custom_parser(x):
	return datetime.strptime('195'+x, '%Y-%m')

# load dataset
series = read_csv('shampoo.csv', header=0, index_col=0, date_parser=custom_parser)

数据集有三年，或 36 个观测值。我们将使用前 24 个用于训练，其余 12 个用作测试集。

下面列出了搜索洗发水销售单变量时间序列预测问题的完整示例网格。

# grid search ets models for monthly shampoo sales
from math import sqrt
from multiprocessing import cpu_count
from joblib import Parallel
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.metrics import mean_squared_error
from pandas import read_csv
from numpy import array

# one-step Holt Winter’s Exponential Smoothing forecast
def exp_smoothing_forecast(history, config):
	t,d,s,p,b,r = config
	# define model
	history = array(history)
	model = ExponentialSmoothing(history, trend=t, damped=d, seasonal=s, seasonal_periods=p)
	# fit model
	model_fit = model.fit(optimized=True, use_boxcox=b, remove_bias=r)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = exp_smoothing_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

# create a set of exponential smoothing configs to try
def exp_smoothing_configs(seasonal=[None]):
	models = list()
	# define config lists
	t_params = ['add', 'mul', None]
	d_params = [True, False]
	s_params = ['add', 'mul', None]
	p_params = seasonal
	b_params = [True, False]
	r_params = [True, False]
	# create config instances
	for t in t_params:
		for d in d_params:
			for s in s_params:
				for p in p_params:
					for b in b_params:
						for r in r_params:
							cfg = [t,d,s,p,b,r]
							models.append(cfg)
	return models

if`_name_`== '__main__':
	# load dataset
	series = read_csv('shampoo.csv', header=0, index_col=0)
	data = series.values
	# data split
	n_test = 12
	# model configs
	cfg_list = exp_smoothing_configs()
	# grid search
	scores = grid_search(data[:,0], cfg_list, n_test)
	print('done')
	# list top 3 configs
	for cfg, error in scores[:3]:
		print(cfg, error)

鉴于存在少量观察，运行该示例很快。

在评估模型时打印模型配置和 RMSE。在运行结束时报告前三个模型配置及其错误。

我们可以看到最好的结果是 RMSE 约为 83.74 销售，具有以下配置：

• 趋势：乘法
• 阻尼：错误
• 季节性：无
• 季节性时期：无
• Box-Cox 变换：错误
• 删除偏差：错误

 > Model[['add', False, None, None, False, True]] 106.431
 > Model[['add', False, None, None, False, False]] 104.874
 > Model[['add', True, None, None, False, False]] 103.069
 > Model[['add', True, None, None, False, True]] 97.918
 > Model[['mul', True, None, None, False, True]] 95.337
 > Model[['mul', True, None, None, False, False]] 102.152
 > Model[['mul', False, None, None, False, True]] 86.406
 > Model[['mul', False, None, None, False, False]] 83.747
 > Model[[None, False, None, None, False, True]] 99.416
 > Model[[None, False, None, None, False, False]] 108.031
done

['mul', False, None, None, False, False] 83.74666940175238
['mul', False, None, None, False, True] 86.40648953786152
['mul', True, None, None, False, True] 95.33737598817238

案例研究 3：季节性

“月平均温度”数据集总结了 1920 至 1939 年华氏诺丁汉城堡的月平均气温，以华氏度为单位。

数据集具有明显的季节性成分，没有明显的趋势。

月平均气温数据集的线图

您可以从 DataMarket 了解有关数据集的更多信息。

直接从这里下载数据集：

• monthly-mean-temp.csv

在当前工作目录中使用文件名“monthly-mean-temp.csv”保存文件。

我们可以使用函数read_csv()将此数据集作为 Pandas 系列加载。

series = read_csv('monthly-mean-temp.csv', header=0, index_col=0)

数据集有 20 年，或 240 个观测值。

我们将数据集修剪为过去五年的数据（60 个观测值），以加快模型评估过程，并使用去年或 12 个观测值来测试集。

# trim dataset to 5 years
data = data[-(5*12):]

季节性成分的周期约为一年，或 12 个观测值。

在准备模型配置时，我们将此作为调用exp_smoothing_configs()函数的季节性时段。

# model configs
cfg_list = exp_smoothing_configs(seasonal=[0, 12])

下面列出了搜索月平均温度时间序列预测问题的完整示例网格。

# grid search ets hyperparameters for monthly mean temp dataset
from math import sqrt
from multiprocessing import cpu_count
from joblib import Parallel
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.metrics import mean_squared_error
from pandas import read_csv
from numpy import array

# one-step Holt Winter’s Exponential Smoothing forecast
def exp_smoothing_forecast(history, config):
	t,d,s,p,b,r = config
	# define model
	history = array(history)
	model = ExponentialSmoothing(history, trend=t, damped=d, seasonal=s, seasonal_periods=p)
	# fit model
	model_fit = model.fit(optimized=True, use_boxcox=b, remove_bias=r)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = exp_smoothing_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

# create a set of exponential smoothing configs to try
def exp_smoothing_configs(seasonal=[None]):
	models = list()
	# define config lists
	t_params = ['add', 'mul', None]
	d_params = [True, False]
	s_params = ['add', 'mul', None]
	p_params = seasonal
	b_params = [True, False]
	r_params = [True, False]
	# create config instances
	for t in t_params:
		for d in d_params:
			for s in s_params:
				for p in p_params:
					for b in b_params:
						for r in r_params:
							cfg = [t,d,s,p,b,r]
							models.append(cfg)
	return models

if`_name_`== '__main__':
	# load dataset
	series = read_csv('monthly-mean-temp.csv', header=0, index_col=0)
	data = series.values
	# trim dataset to 5 years
	data = data[-(5*12):]
	# data split
	n_test = 12
	# model configs
	cfg_list = exp_smoothing_configs(seasonal=[0,12])
	# grid search
	scores = grid_search(data[:,0], cfg_list, n_test)
	print('done')
	# list top 3 configs
	for cfg, error in scores[:3]:
		print(cfg, error)

鉴于大量数据，运行示例相对较慢。

在评估模型时打印模型配置和 RMSE。在运行结束时报告前三个模型配置及其错误。

我们可以看到最好的结果是大约 1.50 度的 RMSE，具有以下配置：

• 趋势：无
• 阻尼：错误
• 季节性：添加剂
• 季节性时期：12
• Box-Cox 变换：错误
• 删除偏差：错误

 > Model[['add', True, 'mul', 12, True, False]] 1.659
 > Model[['add', True, 'mul', 12, True, True]] 1.663
 > Model[['add', True, 'mul', 12, False, True]] 1.603
 > Model[['add', True, 'mul', 12, False, False]] 1.609
 > Model[['mul', False, None, 0, True, True]] 4.920
 > Model[['mul', False, None, 0, True, False]] 4.881
 > Model[['mul', False, None, 0, False, True]] 4.838
 > Model[['mul', False, None, 0, False, False]] 4.813
 > Model[['add', True, 'add', 12, False, True]] 1.568
 > Model[['mul', False, None, 12, True, True]] 4.920
 > Model[['add', True, 'add', 12, False, False]] 1.555
 > Model[['add', True, 'add', 12, True, False]] 1.638
 > Model[['add', True, 'add', 12, True, True]] 1.646
 > Model[['mul', False, None, 12, True, False]] 4.881
 > Model[['mul', False, None, 12, False, True]] 4.838
 > Model[['mul', False, None, 12, False, False]] 4.813
 > Model[['add', True, None, 0, True, True]] 4.654
 > Model[[None, False, 'add', 12, True, True]] 1.508
 > Model[['add', True, None, 0, True, False]] 4.597
 > Model[['add', True, None, 0, False, True]] 4.800
 > Model[[None, False, 'add', 12, True, False]] 1.507
 > Model[['add', True, None, 0, False, False]] 4.760
 > Model[[None, False, 'add', 12, False, True]] 1.502
 > Model[['add', True, None, 12, True, True]] 4.654
 > Model[[None, False, 'add', 12, False, False]] 1.502
 > Model[['add', True, None, 12, True, False]] 4.597
 > Model[[None, False, 'mul', 12, True, True]] 1.507
 > Model[['add', True, None, 12, False, True]] 4.800
 > Model[[None, False, 'mul', 12, True, False]] 1.507
 > Model[['add', True, None, 12, False, False]] 4.760
 > Model[[None, False, 'mul', 12, False, True]] 1.502
 > Model[['add', False, 'add', 12, True, True]] 1.859
 > Model[[None, False, 'mul', 12, False, False]] 1.502
 > Model[[None, False, None, 0, True, True]] 5.188
 > Model[[None, False, None, 0, True, False]] 5.143
 > Model[[None, False, None, 0, False, True]] 5.187
 > Model[[None, False, None, 0, False, False]] 5.143
 > Model[[None, False, None, 12, True, True]] 5.188
 > Model[[None, False, None, 12, True, False]] 5.143
 > Model[[None, False, None, 12, False, True]] 5.187
 > Model[[None, False, None, 12, False, False]] 5.143
 > Model[['add', False, 'add', 12, True, False]] 1.825
 > Model[['add', False, 'add', 12, False, True]] 1.706
 > Model[['add', False, 'add', 12, False, False]] 1.710
 > Model[['add', False, 'mul', 12, True, True]] 1.882
 > Model[['add', False, 'mul', 12, True, False]] 1.739
 > Model[['add', False, 'mul', 12, False, True]] 1.580
 > Model[['add', False, 'mul', 12, False, False]] 1.581
 > Model[['add', False, None, 0, True, True]] 4.980
 > Model[['add', False, None, 0, True, False]] 4.900
 > Model[['add', False, None, 0, False, True]] 5.203
 > Model[['add', False, None, 0, False, False]] 5.151
 > Model[['add', False, None, 12, True, True]] 4.980
 > Model[['add', False, None, 12, True, False]] 4.900
 > Model[['add', False, None, 12, False, True]] 5.203
 > Model[['add', False, None, 12, False, False]] 5.151
 > Model[['mul', True, 'add', 12, True, True]] 19.353
 > Model[['mul', True, 'add', 12, True, False]] 9.807
 > Model[['mul', True, 'add', 12, False, True]] 11.696
 > Model[['mul', True, 'add', 12, False, False]] 2.847
 > Model[['mul', True, None, 0, True, True]] 4.607
 > Model[['mul', True, None, 0, True, False]] 4.570
 > Model[['mul', True, None, 0, False, True]] 4.630
 > Model[['mul', True, None, 0, False, False]] 4.596
 > Model[['mul', True, None, 12, True, True]] 4.607
 > Model[['mul', True, None, 12, True, False]] 4.570
 > Model[['mul', True, None, 12, False, True]] 4.630
 > Model[['mul', True, None, 12, False, False]] 4.593
 > Model[['mul', False, 'add', 12, True, True]] 4.230
 > Model[['mul', False, 'add', 12, True, False]] 4.157
 > Model[['mul', False, 'add', 12, False, True]] 1.538
 > Model[['mul', False, 'add', 12, False, False]] 1.520
done

[None, False, 'add', 12, False, False] 1.5015527325330889
[None, False, 'add', 12, False, True] 1.5015531225114707
[None, False, 'mul', 12, False, False] 1.501561363221282

案例研究 4：趋势和季节性

“月度汽车销售”数据集总结了 1960 年至 1968 年间加拿大魁北克省的月度汽车销量。

数据集具有明显的趋势和季节性成分。

月度汽车销售数据集的线图

您可以从 DataMarket 了解有关数据集的更多信息。

直接从这里下载数据集：

• month-car-sales.csv

在当前工作目录中使用文件名“monthly-car-sales.csv”保存文件。

我们可以使用函数read_csv()将此数据集作为 Pandas 系列加载。

series = read_csv('monthly-car-sales.csv', header=0, index_col=0)

数据集有九年，或 108 个观测值。我们将使用去年或 12 个观测值作为测试集。

季节性成分的期限可能是六个月或 12 个月。在准备模型配置时，我们将尝试将两者作为调用exp_smoothing_configs()函数的季节性时段。

# model configs
cfg_list = exp_smoothing_configs(seasonal=[0,6,12])

下面列出了搜索月度汽车销售时间序列预测问题的完整示例网格。

# grid search ets models for monthly car sales
from math import sqrt
from multiprocessing import cpu_count
from joblib import Parallel
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.metrics import mean_squared_error
from pandas import read_csv
from numpy import array

# one-step Holt Winter’s Exponential Smoothing forecast
def exp_smoothing_forecast(history, config):
	t,d,s,p,b,r = config
	# define model
	history = array(history)
	model = ExponentialSmoothing(history, trend=t, damped=d, seasonal=s, seasonal_periods=p)
	# fit model
	model_fit = model.fit(optimized=True, use_boxcox=b, remove_bias=r)
	# make one step forecast
	yhat = model_fit.predict(len(history), len(history))
	return yhat[0]

# root mean squared error or rmse
def measure_rmse(actual, predicted):
	return sqrt(mean_squared_error(actual, predicted))

# split a univariate dataset into train/test sets
def train_test_split(data, n_test):
	return data[:-n_test], data[-n_test:]

# walk-forward validation for univariate data
def walk_forward_validation(data, n_test, cfg):
	predictions = list()
	# split dataset
	train, test = train_test_split(data, n_test)
	# seed history with training dataset
	history = [x for x in train]
	# step over each time-step in the test set
	for i in range(len(test)):
		# fit model and make forecast for history
		yhat = exp_smoothing_forecast(history, cfg)
		# store forecast in list of predictions
		predictions.append(yhat)
		# add actual observation to history for the next loop
		history.append(test[i])
	# estimate prediction error
	error = measure_rmse(test, predictions)
	return error

# score a model, return None on failure
def score_model(data, n_test, cfg, debug=False):
	result = None
	# convert config to a key
	key = str(cfg)
	# show all warnings and fail on exception if debugging
	if debug:
		result = walk_forward_validation(data, n_test, cfg)
	else:
		# one failure during model validation suggests an unstable config
		try:
			# never show warnings when grid searching, too noisy
			with catch_warnings():
				filterwarnings("ignore")
				result = walk_forward_validation(data, n_test, cfg)
		except:
			error = None
	# check for an interesting result
	if result is not None:
		print(' > Model[%s] %.3f' % (key, result))
	return (key, result)

# grid search configs
def grid_search(data, cfg_list, n_test, parallel=True):
	scores = None
	if parallel:
		# execute configs in parallel
		executor = Parallel(n_jobs=cpu_count(), backend='multiprocessing')
		tasks = (delayed(score_model)(data, n_test, cfg) for cfg in cfg_list)
		scores = executor(tasks)
	else:
		scores = [score_model(data, n_test, cfg) for cfg in cfg_list]
	# remove empty results
	scores = [r for r in scores if r[1] != None]
	# sort configs by error, asc
	scores.sort(key=lambda tup: tup[1])
	return scores

# create a set of exponential smoothing configs to try
def exp_smoothing_configs(seasonal=[None]):
	models = list()
	# define config lists
	t_params = ['add', 'mul', None]
	d_params = [True, False]
	s_params = ['add', 'mul', None]
	p_params = seasonal
	b_params = [True, False]
	r_params = [True, False]
	# create config instances
	for t in t_params:
		for d in d_params:
			for s in s_params:
				for p in p_params:
					for b in b_params:
						for r in r_params:
							cfg = [t,d,s,p,b,r]
							models.append(cfg)
	return models

if`_name_`== '__main__':
	# load dataset
	series = read_csv('monthly-car-sales.csv', header=0, index_col=0)
	data = series.values
	# data split
	n_test = 12
	# model configs
	cfg_list = exp_smoothing_configs(seasonal=[0,6,12])
	# grid search
	scores = grid_search(data[:,0], cfg_list, n_test)
	print('done')
	# list top 3 configs
	for cfg, error in scores[:3]:
		print(cfg, error)

鉴于大量数据，运行示例很慢。

在评估模型时打印模型配置和 RMSE。在运行结束时报告前三个模型配置及其错误。

我们可以看到最好的结果是具有以下配置的约 1,672 销售额的 RMSE：

• 趋势：添加剂
• 阻尼：错误
• 季节性：添加剂
• 季节性时期：12
• Box-Cox 变换：错误
• 删除偏差：是的

这有点令人惊讶，因为我猜想六个月的季节性模型将是首选方法。

 > Model[['add', True, 'add', 6, False, True]] 3240.433
 > Model[['add', True, 'add', 6, False, False]] 3226.384
 > Model[['add', True, 'add', 6, True, False]] 2836.535
 > Model[['add', True, 'add', 6, True, True]] 2784.852
 > Model[['add', True, 'add', 12, False, False]] 1696.173
 > Model[['add', True, 'add', 12, False, True]] 1721.746
 > Model[[None, False, 'add', 6, True, True]] 3204.874
 > Model[['add', True, 'add', 12, True, False]] 2064.937
 > Model[['add', True, 'add', 12, True, True]] 2098.844
 > Model[[None, False, 'add', 6, True, False]] 3190.972
 > Model[[None, False, 'add', 6, False, True]] 3147.623
 > Model[[None, False, 'add', 6, False, False]] 3126.527
 > Model[[None, False, 'add', 12, True, True]] 1834.910
 > Model[[None, False, 'add', 12, True, False]] 1872.081
 > Model[[None, False, 'add', 12, False, True]] 1736.264
 > Model[[None, False, 'add', 12, False, False]] 1807.325
 > Model[[None, False, 'mul', 6, True, True]] 2993.566
 > Model[[None, False, 'mul', 6, True, False]] 2979.123
 > Model[[None, False, 'mul', 6, False, True]] 3025.876
 > Model[[None, False, 'mul', 6, False, False]] 3009.999
 > Model[['add', True, 'mul', 6, True, True]] 2956.728
 > Model[[None, False, 'mul', 12, True, True]] 1972.547
 > Model[[None, False, 'mul', 12, True, False]] 1989.234
 > Model[[None, False, 'mul', 12, False, True]] 1925.010
 > Model[[None, False, 'mul', 12, False, False]] 1941.217
 > Model[[None, False, None, 0, True, True]] 3801.741
 > Model[[None, False, None, 0, True, False]] 3783.966
 > Model[[None, False, None, 0, False, True]] 3801.560
 > Model[[None, False, None, 0, False, False]] 3783.966
 > Model[[None, False, None, 6, True, True]] 3801.741
 > Model[[None, False, None, 6, True, False]] 3783.966
 > Model[[None, False, None, 6, False, True]] 3801.560
 > Model[[None, False, None, 6, False, False]] 3783.966
 > Model[[None, False, None, 12, True, True]] 3801.741
 > Model[[None, False, None, 12, True, False]] 3783.966
 > Model[[None, False, None, 12, False, True]] 3801.560
 > Model[[None, False, None, 12, False, False]] 3783.966
 > Model[['add', True, 'mul', 6, True, False]] 2932.827
 > Model[['mul', True, 'mul', 12, True, True]] 1953.405
 > Model[['add', True, 'mul', 6, False, True]] 2997.259
 > Model[['mul', True, 'mul', 12, True, False]] 1960.242
 > Model[['add', True, 'mul', 6, False, False]] 2979.248
 > Model[['mul', True, 'mul', 12, False, True]] 1907.792
 > Model[['add', True, 'mul', 12, True, True]] 1972.550
 > Model[['add', True, 'mul', 12, True, False]] 1989.236
 > Model[['mul', True, None, 0, True, True]] 3951.024
 > Model[['mul', True, None, 0, True, False]] 3930.394
 > Model[['mul', True, None, 0, False, True]] 3947.281
 > Model[['mul', True, None, 0, False, False]] 3926.082
 > Model[['mul', True, None, 6, True, True]] 3951.026
 > Model[['mul', True, None, 6, True, False]] 3930.389
 > Model[['mul', True, None, 6, False, True]] 3946.654
 > Model[['mul', True, None, 6, False, False]] 3926.026
 > Model[['mul', True, None, 12, True, True]] 3951.027
 > Model[['mul', True, None, 12, True, False]] 3930.368
 > Model[['mul', True, None, 12, False, True]] 3942.037
 > Model[['mul', True, None, 12, False, False]] 3920.756
 > Model[['add', True, 'mul', 12, False, True]] 1750.480
 > Model[['mul', False, 'add', 6, True, False]] 5043.557
 > Model[['mul', False, 'add', 6, False, True]] 7425.711
 > Model[['mul', False, 'add', 6, False, False]] 7448.455
 > Model[['mul', False, 'add', 12, True, True]] 2160.794
 > Model[['mul', False, 'add', 12, True, False]] 2346.478
 > Model[['mul', False, 'add', 12, False, True]] 16303.868
 > Model[['mul', False, 'add', 12, False, False]] 10268.636
 > Model[['mul', False, 'mul', 12, True, True]] 3012.036
 > Model[['mul', False, 'mul', 12, True, False]] 3005.824
 > Model[['add', True, 'mul', 12, False, False]] 1774.636
 > Model[['mul', False, 'mul', 12, False, True]] 14676.476
 > Model[['add', True, None, 0, True, True]] 3935.674
 > Model[['mul', False, 'mul', 12, False, False]] 13988.754
 > Model[['mul', False, None, 0, True, True]] 3804.906
 > Model[['mul', False, None, 0, True, False]] 3805.342
 > Model[['mul', False, None, 0, False, True]] 3778.444
 > Model[['mul', False, None, 0, False, False]] 3798.003
 > Model[['mul', False, None, 6, True, True]] 3804.906
 > Model[['mul', False, None, 6, True, False]] 3805.342
 > Model[['mul', False, None, 6, False, True]] 3778.456
 > Model[['mul', False, None, 6, False, False]] 3798.007
 > Model[['add', True, None, 0, True, False]] 3915.499
 > Model[['mul', False, None, 12, True, True]] 3804.906
 > Model[['mul', False, None, 12, True, False]] 3805.342
 > Model[['mul', False, None, 12, False, True]] 3778.457
 > Model[['mul', False, None, 12, False, False]] 3797.989
 > Model[['add', True, None, 0, False, True]] 3924.442
 > Model[['add', True, None, 0, False, False]] 3905.627
 > Model[['add', True, None, 6, True, True]] 3935.658
 > Model[['add', True, None, 6, True, False]] 3913.420
 > Model[['add', True, None, 6, False, True]] 3924.287
 > Model[['add', True, None, 6, False, False]] 3913.618
 > Model[['add', True, None, 12, True, True]] 3935.673
 > Model[['add', True, None, 12, True, False]] 3913.428
 > Model[['add', True, None, 12, False, True]] 3924.487
 > Model[['add', True, None, 12, False, False]] 3913.529
 > Model[['add', False, 'add', 6, True, True]] 3220.532
 > Model[['add', False, 'add', 6, True, False]] 3199.766
 > Model[['add', False, 'add', 6, False, True]] 3243.478
 > Model[['add', False, 'add', 6, False, False]] 3226.955
 > Model[['add', False, 'add', 12, True, True]] 1833.481
 > Model[['add', False, 'add', 12, True, False]] 1833.511
 > Model[['add', False, 'add', 12, False, True]] 1672.554
 > Model[['add', False, 'add', 12, False, False]] 1680.845
 > Model[['add', False, 'mul', 6, True, True]] 3014.447
 > Model[['add', False, 'mul', 6, True, False]] 3016.207
 > Model[['add', False, 'mul', 6, False, True]] 3025.870
 > Model[['add', False, 'mul', 6, False, False]] 3010.015
 > Model[['add', False, 'mul', 12, True, True]] 1982.087
 > Model[['add', False, 'mul', 12, True, False]] 1981.089
 > Model[['add', False, 'mul', 12, False, True]] 1898.045
 > Model[['add', False, 'mul', 12, False, False]] 1894.397
 > Model[['add', False, None, 0, True, True]] 3815.765
 > Model[['add', False, None, 0, True, False]] 3813.234
 > Model[['add', False, None, 0, False, True]] 3805.649
 > Model[['add', False, None, 0, False, False]] 3809.864
 > Model[['add', False, None, 6, True, True]] 3815.765
 > Model[['add', False, None, 6, True, False]] 3813.234
 > Model[['add', False, None, 6, False, True]] 3805.619
 > Model[['add', False, None, 6, False, False]] 3809.846
 > Model[['add', False, None, 12, True, True]] 3815.765
 > Model[['add', False, None, 12, True, False]] 3813.234
 > Model[['add', False, None, 12, False, True]] 3805.638
 > Model[['add', False, None, 12, False, False]] 3809.837
 > Model[['mul', True, 'add', 6, True, False]] 4099.032
 > Model[['mul', True, 'add', 6, False, True]] 3818.567
 > Model[['mul', True, 'add', 6, False, False]] 3745.142
 > Model[['mul', True, 'add', 12, True, True]] 2203.354
 > Model[['mul', True, 'add', 12, True, False]] 2284.172
 > Model[['mul', True, 'add', 12, False, True]] 2842.605
 > Model[['mul', True, 'add', 12, False, False]] 2086.899
done

['add', False, 'add', 12, False, True] 1672.5539372356582
['add', False, 'add', 12, False, False] 1680.845043013083
['add', True, 'add', 12, False, False] 1696.1734099400082

扩展

本节列出了一些扩展您可能希望探索的教程的想法。

• 数据转换。更新框架以支持可配置的数据转换，例如规范化和标准化。
• 地块预测。更新框架以重新拟合具有最佳配置的模型并预测整个测试数据集，然后将预测与测试集中的实际观察值进行比较。
• 调整历史数量。更新框架以调整用于拟合模型的历史数据量（例如，在 10 年最高温度数据的情况下）。

如果你探索任何这些扩展，我很想知道。

进一步阅读

如果您希望深入了解，本节将提供有关该主题的更多资源。

图书

• 第七章指数平滑，预测：原则和实践，2013。
• 第 6.4 节。时间序列分析简介，工程统计手册，2012。
• 实际时间序列预测与 R ，2016 年。

蜜蜂

• statsmodels.tsa.holtwinters.ExponentialSmoothing API
• statsmodels.tsa.holtwinters.HoltWintersResults API
• Joblib：运行 Python 函数作为管道作业

用品

• 维基百科上的指数平滑

摘要

在本教程中，您了解了如何开发一个框架，用于网格搜索所有指数平滑模型超参数，以进行单变量时间序列预测。

具体来说，你学到了：

• 如何使用前向验证从零开始开发网格搜索 ETS 模型的框架。
• 如何为出生日常时间序列数据网格搜索 ETS 模型超参数。
• 如何为洗发水销售，汽车销售和温度的月度时间序列数据网格搜索 ETS 模型超参数。

你有任何问题吗？ 在下面的评论中提出您的问题，我会尽力回答。

一个标准的人类活动识别问题的温和介绍

原文： machinelearningmastery.com/how-to-load-and-explore-a-standard-human-activity-recognition-problem/

人类活动识别是将由专用线束或智能电话记录的加速度计数据序列分类为已知的明确定义的运动的问题。

鉴于每秒产生大量观测结果，观测的时间性质以及缺乏将加速度计数据与已知运动联系起来的明确方法，这是一个具有挑战性的问题。

该问题的经典方法涉及来自基于固定大小窗口的时间序列数据的手工制作特征和训练机器学习模型，例如决策树的集合。困难在于此功能工程需要该领域的深厚专业知识。

最近，诸如循环神经网络和一维卷积神经网络或 CNN 的深度学习方法已经被证明在很少或没有数据特征工程的情况下提供具有挑战性的活动识别任务的最新结果。

在本教程中，您将发现用于时间序列分类的标准人类活动识别数据集以及如何在建模之前探索数据集。

完成本教程后，您将了解：

• 如何加载和准备人类活动识别时间序列分类数据。
• 如何探索和可视化时间序列分类数据，以生成建模的想法。
• 一套用于构建问题，准备数据，建模和评估人类活动识别模型的方法。

让我们开始吧。

标准人类活动识别问题的温和介绍 tgraham 的照片，保留一些权利。

教程概述

本教程分为 8 个部分;他们是：

1. 人类活动识别
2. 问题描述
3. 加载数据集
4. 绘制单个主题的跟踪
5. 绘制总活动持续时间
6. 绘制每个主题的痕迹
7. 绘制痕量观测的直方图
8. 建模问题的方法

人类活动识别

人类活动识别或简称 HAR，是基于使用传感器的移动痕迹来预测人正在做什么的问题。

运动通常是正常的室内活动，例如站立，坐着，跳跃和上楼梯。传感器通常位于主体上，例如智能手机或背心，并且经常以三维（x，y，z）记录加速度计数据。

这个想法是，一旦主体的活动被识别和知道，智能计算机系统就可以提供帮助。

这是一个具有挑战性的问题，因为没有明确的分析方法将传感器数据与一般方式的特定动作联系起来。由于收集了大量的传感器数据（例如，每秒数十或数百次观察），并且在开发预测模型时从这些数据中经典使用手工制作的特征和启发式，因此在技术上具有挑战性。

最近，深度学习方法已经在 HAR 问题上取得了成功，因为它们能够自动学习更高阶的特征。

基于传感器的活动识别从大量低水平传感器读数中寻找关于人类活动的深刻的高级知识。传统的模式识别方法在过去几年中取得了巨大的进步。然而，这些方法通常严重依赖于启发式手工特征提取，这可能会妨碍它们的泛化表现。 [...]最近，深度学习的最新进展使得可以执行自动高级特征提取，从而在许多领域实现了有希望的表现。

• 基于传感器的活动识别深度学习：调查

问题描述

收集数据集“来自单个胸部安装的加速度计数据集的 _ 活动识别”并由 Casale，Pujol 等人提供。来自西班牙巴塞罗那大学。_

它可以从 UCI 机器学习库免费获得：

• 来自单个胸部加速计数据集，UCI 机器学习库的活动识别。

在 2011 年论文“使用可穿戴设备的加速度计数据人类活动识别”中描述并使用该数据集作为序列分类模型的基础。

数据集由来自 15 个不同科目的未校准加速度计数据组成，每个科目执行 7 项活动。每个受试者佩戴定制的胸部加速度计，并以 52Hz（每秒 52 次观察）收集数据。

每个主题佩戴的定制胸部系统的照片。 取自“使用可穿戴设备从加速度计数据中识别人体活动”。

每个受试者以连续的顺序进行和记录七项活动。

所开展的具体活动包括：

• 1：在计算机工作（workingPC）。
• 2：站起来，走路，上下楼梯。
• 3：站立（站立）。
• 4：散步（醒来）。
• 5：上/下楼梯（楼梯）。
• 6：与某人散步和交谈。
• 7：站立时说话（说话）。

该论文使用了数据的一个子集，特别是 14 个主题和 5 个活动。目前尚不清楚为什么没有使用其他 2 项活动（2 和 6）。

数据来自 14 名测试者，3 名女性和 11 名年龄在 27 岁至 35 岁之间的男性。[...]所收集的数据由上下楼梯步行 33 分钟，步行 82 分钟，说话 115 分钟组成。保持站立的分钟数和 86 分钟的电脑工作时间。

• 使用可穿戴设备从加速计数据识别人类活动，2011。

本文的重点是从数据开发手工制作的功能，以及机器学习算法的开发。

数据被分成一秒钟的观察窗口（52），窗口之间有 50％的重叠。

我们使用 52 个样本的窗口从数据中提取特征，对应于 1 秒的加速度计数据，窗口之间有 50％的重叠。从每个窗口，我们建议提取以下特征：窗口中加速度积分的均方根值，以及 Minmax 和的平均值。 [...]尽管如此，为了完成这组特征，我们添加了已被证明对人类活动识别有用的特征，如：平均值，标准偏差，偏度，峰度，加速度计轴的每对成对之间的相关性（不包括幅度） ，七级小波分解系数的能量。通过这种方式，我们获得了 319 维特征向量。

• 使用可穿戴设备从加速计数据识别人类活动，2011。

使用 5 倍交叉验证拟合和评估一套机器学习模型，并获得 94％的准确度。

随机森林的混淆矩阵在数据集上进行评估。 取自“使用可穿戴设备从加速度计数据中识别人体活动”。

加载数据集

数据集可以直接从 UCI 机器学习库下载。

• 下载数据集

将数据集解压缩到名为“HAR”的新目录中。

该目录包含 CSV 文件列表，每个主题一个（1-15）和自述文件。

每个文件包含 5 列，行号，x，y 和 z 加速度计读数以及 0 到 7 的类号，其中类 0 表示“无活动”，类 1-7 对应于上一节中列出的活动。

例如，下面是文件“1.csv”的前 5 行：

0,1502,2215,2153,1
1,1667,2072,2047,1
2,1611,1957,1906,1
3,1601,1939,1831,1
4,1643,1965,1879,1
...

首先，我们可以将每个文件作为单个 NumPy 数组加载并删除第一列。

下面名为load_dataset()的函数将加载 HAR 目录中的所有 CSV 文件，删除第一列并返回 15 个 NumPy 数组的列表。

# load sequence for each subject, returns a list of numpy arrays
def load_dataset(prefix=''):
	subjects = list()
	directory = prefix + 'HAR/'
	for name in listdir(directory):
		filename = directory + '/' + name
		if not filename.endswith('.csv'):
			continue
		df = read_csv(filename, header=None)
		# drop row number
		values = df.values[:, 1:]
		subjects.append(values)
	return subjects

下面列出了完整的示例。

# load dataset
from os import listdir
from pandas import read_csv

# load sequence for each subject, returns a list of numpy arrays
def load_dataset(prefix=''):
	subjects = list()
	directory = prefix + 'HAR/'
	for name in listdir(directory):
		filename = directory + '/' + name
		if not filename.endswith('.csv'):
			continue
		df = read_csv(filename, header=None)
		# drop row number
		values = df.values[:, 1:]
		subjects.append(values)
	return subjects

# load
subjects = load_dataset()
print('Loaded %d subjects' % len(subjects))

运行该示例将加载所有数据并打印已加载主题的数量。

Loaded 15 subjects

注意，目录中的文件是按文件顺序导航的，这可能与主题顺序不同，例如10.csv以文件顺序出现在2.csv之前。我相信这在本教程中无关紧要。

现在我们知道了如何加载数据，我们可以通过一些可视化来探索它。

绘制单个主题的跟踪

良好的第一个可视化是绘制单个主题的数据。

我们可以为给定主题的每个变量创建一个图形，包括 x，y 和 z 加速度计数据，以及相关的类类值。

下面的函数plot_subject()将绘制给定主题的数据。

# plot the x, y, z acceleration and activities for a single subject
def plot_subject(subject):
	pyplot.figure()
	# create a plot for each column
	for col in range(subject.shape[1]):
		pyplot.subplot(subject.shape[1], 1, col+1)
		pyplot.plot(subject[:,col])
	pyplot.show()

我们可以将它与上一节中的数据加载结合起来，并绘制第一个加载主题的数据。

# plot a subject
from os import listdir
from pandas import read_csv
from matplotlib import pyplot

# load sequence for each subject, returns a list of numpy arrays
def load_dataset(prefix=''):
	subjects = list()
	directory = prefix + 'HAR/'
	for name in listdir(directory):
		filename = directory + '/' + name
		if not filename.endswith('.csv'):
			continue
		df = read_csv(filename, header=None)
		# drop row number
		values = df.values[:, 1:]
		subjects.append(values)
	return subjects

# plot the x, y, z acceleration and activities for a single subject
def plot_subject(subject):
	pyplot.figure()
	# create a plot for each column
	for col in range(subject.shape[1]):
		pyplot.subplot(subject.shape[1], 1, col+1)
		pyplot.plot(subject[:,col])
	pyplot.show()

# load
subjects = load_dataset()
print('Loaded %d subjects' % len(subjects))
# plot activities for a single subject
plot_subject(subjects[0])

运行该示例为第一个加载主题的每个变量创建一个线图。

我们可以在序列的开头看到一些非常大的运动，这可能是一个可以被删除的异常或异常行为。

我们还可以看到主题多次执行某些操作。例如，仔细查看类变量的图（底部图）表明受试者按以下顺序执行活动：1,2,0,3,0,4,3,5,3,6,7。活动 3 进行了两次。

第一个加载主题的 x，y，z 和类的线图。

我们可以重新运行此代码并绘制第二个加载的主题（可能是10.csv）。

...
# plot activities for a single subject
plot_subject(subjects[1])

运行该示例会创建一个类似的图。

我们可以看到更多细节，这表明在上一个情节开头看到的大异常值可能是从该主题的痕迹中清除了值。

我们看到类似的活动序列，活动 3 再次发生两次。

我们还可以看到，某些活动的执行时间比其他活动长得多。这可能会影响模型区分活动的能力，例如：两个受试者的活动 3（站立）相对于所进行的其他活动具有非常少的数据。

第二个加载主题的 x，y，z 和类的线图。

绘制总活动持续时间

上一节提出了一个问题，即我们对所有主题的每项活动进行了多长时间或多少次观察。

如果一项活动的数据多于另一项活动，这可能很重要，这表明不太好的活动可能难以建模。

我们可以通过按活动和主题对所有观察进行分组来研究这一点，并绘制分布图。这将了解每个主题在跟踪过程中花费多长时间执行每项活动。

首先，我们可以为每个主题分组活动。

我们可以通过为每个主题创建字典并按活动存储所有跟踪数据来完成此操作。下面的group_by_activity()功能将为每个主题执行此分组。

# returns a list of dict, where each dict has one sequence per activity
def group_by_activity(subjects, activities):
	grouped = [{a:s[s[:,-1]==a] for a in activities} for s in subjects]
	return grouped

接下来，我们可以计算每个主题的每个活动的总持续时间。

我们知道加速度计数据是以 52Hz 记录的，因此我们可以将每个活动的每个跟踪的长度除以 52，以便以秒为单位总结持续时间。

以下名为plot_durations()的函数将计算每个主题的每个活动的持续时间，并将结果绘制为箱线图。盒状和须状图是总结每个活动的 15 个持续时间的有用方式，因为它描述了持续时间的扩展而不假设分布。

# calculate total duration in sec for each activity per subject and plot
def plot_durations(grouped, activities):
	# calculate the lengths for each activity for each subject
	freq = 52
	durations = [[len(s[a])/freq for s in grouped] for a in activities]
	pyplot.boxplot(durations, labels=activities)
	pyplot.show()

下面列出了绘制活动持续时间分布的完整示例。

# durations by activity
from os import listdir
from pandas import read_csv
from matplotlib import pyplot

# load sequence for each subject, returns a list of numpy arrays
def load_dataset(prefix=''):
	subjects = list()
	directory = prefix + 'HAR/'
	for name in listdir(directory):
		filename = directory + '/' + name
		if not filename.endswith('.csv'):
			continue
		df = read_csv(filename, header=None)
		# drop row number
		values = df.values[:, 1:]
		subjects.append(values)
	return subjects

# returns a list of dict, where each dict has one sequence per activity
def group_by_activity(subjects, activities):
	grouped = [{a:s[s[:,-1]==a] for a in activities} for s in subjects]
	return grouped

# calculate total duration in sec for each activity per subject and plot
def plot_durations(grouped, activities):
	# calculate the lengths for each activity for each subject
	freq = 52
	durations = [[len(s[a])/freq for s in grouped] for a in activities]
	pyplot.boxplot(durations, labels=activities)
	pyplot.show()

# load
subjects = load_dataset()
print('Loaded %d subjects' % len(subjects))
# group traces by activity for each subject
activities = [i for i in range(0,8)]
grouped = group_by_activity(subjects, activities)
# plot durations
plot_durations(grouped, activities)

运行该示例绘制了每个主题的活动持续时间的分布。

我们可以看到，对于活动 0（无活动），2（站立，行走和上下楼梯），5（上/下楼梯）和 6（步行和说话）的观察相对较少。

这可能表明为什么活动 2 和 6 被排除在原始论文的实验之外。

我们还可以看到每个主题在活动 1（站立，行走和上下楼梯）和活动 7（站立时说话）上花费了大量时间。这些活动可能过多。准备模型数据可能有益于对这些活动进行欠采样或对其他活动进行采样。

每个受试者活动持续时间分布的箱线图

绘制每个主题的痕迹

接下来，查看每个主题的跟踪数据可能会很有趣。

一种方法是在单个图形上绘制单个主题的所有迹线，然后垂直排列所有图形。这将允许跨主题和主题内的痕迹进行比较。

以下名为plot_subjects()的函数将在单独的图上绘制 15 个主题中每个主题的加速度计数据。每个 x，y 和 z 数据的迹线分别绘制为橙色，绿色和蓝色。

# plot the x, y, z acceleration for each subject
def plot_subjects(subjects):
	pyplot.figure()
	# create a plot for each subject
	for i in range(len(subjects)):
		pyplot.subplot(len(subjects), 1, i+1)
		# plot each of x, y and z
		for j in range(subjects[i].shape[1]-1):
			pyplot.plot(subjects[i][:,j])
	pyplot.show()

下面列出了完整的示例。

# plot accelerometer data for all subjects
from os import listdir
from pandas import read_csv
from matplotlib import pyplot

# load sequence for each subject, returns a list of numpy arrays
def load_dataset(prefix=''):
	subjects = list()
	directory = prefix + 'HAR/'
	for name in listdir(directory):
		filename = directory + '/' + name
		if not filename.endswith('.csv'):
			continue
		df = read_csv(filename, header=None)
		# drop row number
		values = df.values[:, 1:]
		subjects.append(values)
	return subjects

# plot the x, y, z acceleration for each subject
def plot_subjects(subjects):
	pyplot.figure()
	# create a plot for each subject
	for i in range(len(subjects)):
		pyplot.subplot(len(subjects), 1, i+1)
		# plot each of x, y and z
		for j in range(subjects[i].shape[1]-1):
			pyplot.plot(subjects[i][:,j])
	pyplot.show()

# load
subjects = load_dataset()
print('Loaded %d subjects' % len(subjects))
# plot trace data for each subject
plot_subjects(subjects)

运行该示例将创建一个包含 15 个图的图形。

我们正在研究一般而非具体的趋势。

• 我们可以看到很多橙色和绿色以及非常小的蓝色，这表明 z 数据在建模这个问题时可能不太重要。
• 我们可以看到跟踪数据在 x 和 y 跟踪的相同时间发生相同的一般变化，这表明可能只需要一个数据轴来拟合预测模型。
• 我们可以看到每个受试者在序列开始时（前 60 秒）在迹线中具有相同的大尖峰，可能与实验的启动有关。
• 我们可以在跟踪数据中看到跨主题的类似结构，尽管某些迹线看起来比其他迹线更柔和，例如比较第一和第二个图上的幅度。

每个受试者可具有不同的完整迹线长度，因此通过 x 轴的直接比较可能是不合理的（例如，同时执行类似的活动）。无论如何，我们并不真正关心这个问题。

所有 15 名受试者的加速度计跟踪数据的线图。

痕迹似乎具有相同的一般比例，但受试者之间的幅度差异表明，每个受试者重新缩放数据可能比跨受试者缩放更有意义。

这对于训练数据可能是有意义的，但对于缩放测试对象的数据可能在方法上具有挑战性。它需要或假设在预测活动之前可以获得整个跟踪。这对于模型的离线使用很好，但不能在线使用模型。它还表明，使用预先校准的跟踪数据（例如以固定规模进入的数据）可以更容易地在线使用模型。

绘制痕量观测的直方图

上一节中关于跨不同主题的显着不同尺度可能性的观点可能会给这个数据集的建模带来挑战。

我们可以通过绘制加速度计数据的每个轴的观测分布的直方图来探索这一点。

与上一节一样，我们可以为每个主题创建一个绘图，然后将所有主题的绘图与相同的 x 轴垂直对齐，以帮助发现展开的明显差异。

更新的plot_subjects()函数用于绘制直方图而不是线图，如下所示。hist()函数用于为加速度计数据的每个轴创建直方图，并且使用大量箱（100）来帮助展开图中的数据。子图也都共享相同的 x 轴以帮助进行比较。

# plot the x, y, z acceleration for each subject
def plot_subjects(subjects):
	pyplot.figure()
	# create a plot for each subject
	xaxis = None
	for i in range(len(subjects)):
		ax = pyplot.subplot(len(subjects), 1, i+1, sharex=xaxis)
		if i == 0:
			xaxis = ax
		# plot a histogram of x data
		for j in range(subjects[i].shape[1]-1):
			pyplot.hist(subjects[i][:,j], bins=100)
	pyplot.show()

下面列出了完整的示例

# plot histograms of trace data for all subjects
from os import listdir
from pandas import read_csv
from matplotlib import pyplot

# load sequence for each subject, returns a list of numpy arrays
def load_dataset(prefix=''):
	subjects = list()
	directory = prefix + 'HAR/'
	for name in listdir(directory):
		filename = directory + '/' + name
		if not filename.endswith('.csv'):
			continue
		df = read_csv(filename, header=None)
		# drop row number
		values = df.values[:, 1:]
		subjects.append(values)
	return subjects

# plot the x, y, z acceleration for each subject
def plot_subjects(subjects):
	pyplot.figure()
	# create a plot for each subject
	xaxis = None
	for i in range(len(subjects)):
		ax = pyplot.subplot(len(subjects), 1, i+1, sharex=xaxis)
		if i == 0:
			xaxis = ax
		# plot a histogram of x data
		for j in range(subjects[i].shape[1]-1):
			pyplot.hist(subjects[i][:,j], bins=100)
	pyplot.show()

# load
subjects = load_dataset()
print('Loaded %d subjects' % len(subjects))
# plot trace data for each subject
plot_subjects(subjects)

运行该示例将创建一个包含 15 个图的单个图形，每个图形对应一个图形，以及每个图表的 3 个加速度计数据的 3 个直方图。

蓝色，橙色和绿色三种颜色代表 x，y 和 z 轴。

该图表明加速度计的每个轴的分布是高斯分布或者非常接近高斯分布。这可以帮助沿着加速度计数据的每个轴进行简单的离群值检测和移除。

该图确实有助于显示主题内的分布与主题之间的分布差异。

在每个主题内，共同模式是 x（蓝色）和 z（绿色）一起分组到左边，y 数据（橙色）分开到右边。 y 的分布通常更尖锐，因为 x 和 z 的分布更平坦。

在整个主题中，我们可以看到一般的聚类值约为 2,000（无论单位是多少），尽管有很多差异。这种显着的分布差异确实表明在执行任何跨主题建模之前，需要至少标准化（转换为零均值和单位方差）每个轴和每个主体的数据。

每个受试者的加速度计数据的直方图

建模问题的方法

在本节中，我们将基于对数据集的上述探索，探索针对该问题的数据准备和建模的一些想法和方法。

这些可能有助于特别是对该数据集建模，但也有助于人类活动识别，甚至是一般的时间序列分类问题。

问题框架

尽管所有方法都围绕时间序列分类的思想，但有许多方法可以将数据构建为预测性建模问题。

要考虑的两种主要方法是：

• 每个受试者：每个受试者的痕量数据的模型活动分类。
• 交叉主题：跨主题的跟踪数据的模型活动分类。

后者，交叉主题，是更理想的，但如果目标是深刻理解给定的主题，例如前者也可能是有趣的。家中的个性化模型。

在建模过程中构建数据的两种主要方法包括：

• 分段活动：跟踪数据可以按活动预先分段，并且可以针对每个活动对整个跟踪或其特征进行训练的模型。
• 滑动窗口：每个主体的连续轨迹被分成滑动窗口，有或没有重叠，窗口的每个活动的模式被视为要预测的活动。

就模型的实际使用而言，前一种方法可能不太现实，但可能是更容易建模的问题。后者是原始论文中使用的问题的框架，其中制备了具有 50％重叠的 1 秒窗口。

我没有看到问题框架中的重叠的好处，除了加倍训练数据集的大小，这可能有益于深度神经网络。事实上，我预计它可能会导致过度模型。

数据准备

数据表明在建模过程中可能有用的一些准备方案：

• 将加速度计观测值下采样到几分之一秒可能是有帮助的，例如， 1 / 4,1 / 2,1,2 秒。
• 截断原始数据的前 60 秒可能是谨慎的，因为它似乎与实验的启动有关，并且所有主体当时正在执行活动 1（在计算机上工作）。
• 使用简单的离群值检测和去除方法（例如，每个轴的平均值的标准差的 3 到 4 倍的值）可能是有用的。
• 也许删除具有相对较少观察的活动将是明智的，或者导致对预测方法（例如，活动 0,2 和 6）的更公平的评估。
• 也许通过对代表性不足的活动进行过度采样或对训练数据集中过度代表的活动进行采样调整来重新平衡活动可能有助于建模。
• 尝试不同的窗口尺寸将是有趣的（例如 1,5,10,30 秒），尤其是在对观察的下采样的确证中。
• 对于任何跨主题模型，几乎肯定需要标准化每个主题的数据。在每个受试者标准化后对受试者的数据进行标准化也可能是有用的。

如前一节所述，每个主题的数据标准化确实引入了方法问题，并且无论如何都可以使用，因为需要来自原始硬件系统的校准观察。

问题建模

我建议使用神经网络来探索这个问题。

与使用特征工程和特定于域的手工制作特征的论文中使用的方法不同，直接对原始数据进行建模（下采样或其他方式）将是有用且通用的。

首先，我建议使用强大的方法（如随机森林或梯度增强机器）发现表现基线。然后探索特别适合时间序列分类问题的神经网络方法。

可能适合的两种类型的神经网络架构是：

• 卷积神经网络或 CNN。
• 循环神经网络或 RNN，特别是长短期记忆或 LSTM。

第三个是两者的混合：

• CNN LSTMs。

CNN 能够从输入序列中提取特征，例如输入加速度计数据的窗口。诸如 LSTM 之类的 RNN 能够直接从长序列的输入数据中学习，并学习数据中的长期关系。

我希望序列数据中几乎没有因果关系，除了每个主题看起来他们正在执行相同的人为行动序列，我们不想学习。朴素地，这可能表明 CNN 更适合于在给定一系列观察到的加速度计数据的情况下预测活动。

一维 CNN 已广泛用于此类问题，其中一个通道用于加速度计数据的每个轴。一个很好的简单起点是直接在序列数据的窗口上拟合 CNN 模型。这是 2014 年题为“使用移动传感器进行人类活动识别的卷积神经网络”的论文中描述的方法，并且从下面的图中可以看出更清楚。

用于人类活动识别的 1D CNN 架构的示例 取自“使用移动传感器的用于人类活动识别的卷积神经网络”。

CNN LSTM 可用于 CNN 学习观察子序列的表示，然后 LSTM 学习这些子序列。

例如，CNN 可以提取一秒钟的加速度计数据，然后可以重复 30 秒，以向 LSTM 提供 30 个 CNN 解释数据的时间步长。

我希望这三个方法对这个问题和类似的问题都很有意思。

模型评估

我不认为窗口的重叠是有用的，并且实际上如果在交叉验证期间跟踪数据的部分在训练和测试数据集中都可用，则实际上可能导致轻微的误导性结果。然而，它会增加训练数据的数量。

我认为重要的是将数据暴露给模型，同时保持观察的时间顺序。对于给定主题，来自单个活动的多个窗口可能看起来相似，并且随机改组和分离窗口以训练测试数据可能导致误导结果。

我不建议在原始论文中使用随机改组的 k-fold 交叉验证。我希望这会带来乐观的结果，每个 15 个主题的痕迹中有一秒钟的数据窗口混合在一起进行训练和测试。

也许对这些数据中的模型进行公平评估将是按主题使用留一法交叉验证或 LOOCV。这是模型适合前 14 个主题并在第 15 个主题的所有窗口上进行评估的地方。重复该过程，其中每个受试者有机会被用作保持测试数据集。

按主题分割数据集避免了在模型评估期间与各个窗口的时间排序相关的任何问题，因为所有窗口都将保证新的/看不见的数据。

如果你探索这些建模思想中的任何一个，我很想知道。

进一步阅读

如果您希望深入了解，本节将提供有关该主题的更多资源。

文件

• 使用加速度计进行活动识别的综合研究，2018
• 基于传感器的活动识别深度学习：调查，2017。
• 使用可穿戴设备从加速度计数据识别人类活动，2011。
• 用于使用移动传感器识别人类活动的卷积神经网络，2014。

API

• sklearn.preprocessing.StandardScaler API
• sklearn.preprocessing.MinMaxScaler API
• keras.utils.to_categorical API
• Keras 卷积层

用品

• 活动识别，维基百科
• 来自单个胸部加速计数据集，UCI 机器学习库的活动识别。

摘要

在本教程中，您发现了一个用于时间序列分类的标准人类活动识别数据集。

具体来说，你学到了：

• 如何加载和准备人类活动识别时间序列分类数据。
• 如何探索和可视化时间序列分类数据，以生成建模的想法。
• 一套用于构建问题，准备数据，建模和评估人类活动识别模型的方法。

你有任何问题吗？ 在下面的评论中提出您的问题，我会尽力回答。

如何加载和探索家庭用电数据

原文： machinelearningmastery.com/how-to-load-and-explore-household-electricity-usage-data/

鉴于智能电表的兴起以及太阳能电池板等发电技术的广泛采用，可提供大量的用电数据。

该数据代表功率相关变量的多变量时间序列，而这些变量又可用于建模甚至预测未来的电力消耗。

在本教程中，您将发现用于多步时间序列预测的家庭功耗数据集，以及如何使用探索性分析更好地理解原始数据。

完成本教程后，您将了解：

• 家庭用电量数据集，描述四年内单个房屋的用电量。
• 如何使用一系列线图来探索和理解数据集，用于数据分布的系列数据和直方图。
• 如何使用对问题的新理解来考虑预测问题的不同框架，可以准备数据的方式以及可以使用的建模方法。

让我们开始吧。

如何加载和探索家庭用电数据 Sheila Sund 的照片，保留一些权利。

教程概述

本教程分为五个部分;他们是：

1. 功耗数据集
2. 加载数据集
3. 随着时间的推移观察模式
4. 时间序列数据分布
5. 关于建模的想法

家庭用电量数据集

家庭用电量数据集是一个多变量时间序列数据集，描述了四年内单个家庭的用电量。

该数据是在 2006 年 12 月至 2010 年 11 月之间收集的，并且每分钟收集家庭内的能耗观察结果。

它是一个多变量系列，由七个变量组成（除日期和时间外）;他们是：

• global_active_power ：家庭消耗的总有功功率（千瓦）。
• global_reactive_power ：家庭消耗的总无功功率（千瓦）。
• 电压：平均电压（伏特）。
• global_intensity ：平均电流强度（安培）。
• sub_metering_1 ：厨房的有功电能（瓦特小时的有功电能）。
• sub_metering_2 ：用于洗衣的有功能量（瓦特小时的有功电能）。
• sub_metering_3 ：气候控制系统的有功电能（瓦特小时的有功电能）。

有功和无功电能参考交流电的技术细节。

一般而言，有功能量是家庭消耗的实际功率，而无功能量是线路中未使用的功率。

我们可以看到，数据集通过房屋中的主电路，特别是厨房，洗衣房和气候控制，提供有功功率以及有功功率的某种划分。这些不是家庭中的所有电路。

通过首先将有功能量转换为瓦特小时，然后以瓦时为单位减去其他亚计量有功能量，可以从有功能量计算剩余瓦特小时，如下所示：

sub_metering_remainder = (global_active_power * 1000 / 60) - (sub_metering_1 + sub_metering_2 + sub_metering_3)

数据集似乎是在没有开创性参考文件的情况下提供的。

尽管如此，该数据集已成为评估多步预测的时间序列预测和机器学习方法的标准，特别是用于预测有功功率。此外，尚不清楚数据集中的其他特征是否可以使模型在预测有功功率方面受益。

加载数据集

数据集可以从 UCI 机器学习库下载为单个 20 兆字节的.zip 文件：

• household_power_consumption.zip

下载数据集并将其解压缩到当前工作目录中。您现在将拥有大约 127 兆字节的文件“household_power_consumption.txt”并包含所有观察结果

检查数据文件。

以下是原始数据文件中的前五行数据（和标题）。

Date;Time;Global_active_power;Global_reactive_power;Voltage;Global_intensity;Sub_metering_1;Sub_metering_2;Sub_metering_3
16/12/2006;17:24:00;4.216;0.418;234.840;18.400;0.000;1.000;17.000
16/12/2006;17:25:00;5.360;0.436;233.630;23.000;0.000;1.000;16.000
16/12/2006;17:26:00;5.374;0.498;233.290;23.000;0.000;2.000;17.000
16/12/2006;17:27:00;5.388;0.502;233.740;23.000;0.000;1.000;17.000
16/12/2006;17:28:00;3.666;0.528;235.680;15.800;0.000;1.000;17.000
...

我们可以看到数据列用分号分隔（'; '）。

据报道，该数据在该时间段内每天有一行。

数据确实缺少值;例如，我们可以在 28/4/2007 左右看到 2-3 天的缺失数据。

...
28/4/2007;00:20:00;0.492;0.208;236.240;2.200;0.000;0.000;0.000
28/4/2007;00:21:00;?;?;?;?;?;?;
28/4/2007;00:22:00;?;?;?;?;?;?;
28/4/2007;00:23:00;?;?;?;?;?;?;
28/4/2007;00:24:00;?;?;?;?;?;?;
...

我们可以通过将数据文件作为 Pandas DataFrame 加载并总结加载的数据来启动。

我们可以使用 read_csv（）函数来加载数据。

使用此功能很容易加载数据，但正确加载它有点棘手。

具体来说，我们需要做一些自定义的事情：

• 将列之间的单独值指定为分号（sep =';'）
• 指定第 0 行具有列的名称（header = 0）
• 指定我们有大量的 RAM 来避免警告我们将数据作为对象数组而不是数组加载，因为缺少数据的'？'值（low_memory = False）。
• 指定 Pandas 在解析日期时尝试推断日期时间格式是可以的，这样会更快（infer_datetime_format = True）
• 指定我们要将日期和时间列一起解析为名为“datetime”的新列（parse_dates = {'datetime'：[0,1]}）
• 指定我们希望新的“datetime”列成为 DataFrame 的索引（index_col = ['datetime']）。

将所有这些放在一起，我们现在可以加载数据并汇总加载的形状和前几行。

# load all data
dataset = read_csv('household_power_consumption.txt', sep=';', header=0, low_memory=False, infer_datetime_format=True, parse_dates={'datetime':[0,1]}, index_col=['datetime'])
# summarize
print(dataset.shape)
print(dataset.head())

接下来，我们可以使用带有 NaN 值的“？”字符标记所有缺失值，这是一个浮点数。

这将允许我们将数据作为一个浮点值数组而不是混合类型来处理，效率较低。

# mark all missing values
dataset.replace('?', nan, inplace=True)

现在，我们可以使用上一节中的计算创建一个包含剩余子计量的新列。

# add a column for for the remainder of sub metering
values = dataset.values.astype('float32')
dataset['sub_metering_4'] = (values[:,0] * 1000 / 60) - (values[:,4] + values[:,5] + values[:,6])

我们现在可以将清理后的数据集版本保存到新文件中;在这种情况下，我们只需将文件扩展名更改为.csv，并将数据集保存为“household_power_consumption.csv”。

# save updated dataset
dataset.to_csv('household_power_consumption.csv')

为了确认我们没有弄乱，我们可以重新加载数据集并汇总前五行。

# load the new file
dataset = read_csv('household_power_consumption.csv', header=None)
print(dataset.head())

将所有这些结合在一起，下面列出了加载，清理和保存数据集的完整示例。

# load and clean-up data
from numpy import nan
from pandas import read_csv
# load all data
dataset = read_csv('household_power_consumption.txt', sep=';', header=0, low_memory=False, infer_datetime_format=True, parse_dates={'datetime':[0,1]}, index_col=['datetime'])
# summarize
print(dataset.shape)
print(dataset.head())
# mark all missing values
dataset.replace('?', nan, inplace=True)
# add a column for for the remainder of sub metering
values = dataset.values.astype('float32')
dataset['sub_metering_4'] = (values[:,0] * 1000 / 60) - (values[:,4] + values[:,5] + values[:,6])
# save updated dataset
dataset.to_csv('household_power_consumption.csv')
# load the new dataset and summarize
dataset = read_csv('household_power_consumption.csv', header=0, infer_datetime_format=True, parse_dates=['datetime'], index_col=['datetime'])
print(dataset.head())

运行该示例首先加载原始数据并汇总加载数据的形状和前五行。

(2075259, 7)

                    Global_active_power      ...       Sub_metering_3
datetime                                     ...
2006-12-16 17:24:00               4.216      ...                 17.0
2006-12-16 17:25:00               5.360      ...                 16.0
2006-12-16 17:26:00               5.374      ...                 17.0
2006-12-16 17:27:00               5.388      ...                 17.0
2006-12-16 17:28:00               3.666      ...                 17.0

然后清理数据集并将其保存到新文件中。

我们加载这个新文件并再次打印前五行，显示删除日期和时间列以及添加新的子计量列。

                     Global_active_power       ...        sub_metering_4
datetime                                       ...
2006-12-16 17:24:00                4.216       ...             52.266670
2006-12-16 17:25:00                5.360       ...             72.333336
2006-12-16 17:26:00                5.374       ...             70.566666
2006-12-16 17:27:00                5.388       ...             71.800000
2006-12-16 17:28:00                3.666       ...             43.100000

我们可以查看新的'household_power_consumption.csv'文件并检查缺失的观察结果是否用空列标记，大熊猫将正确读作 NaN，例如第 190,499 行：

...
2007-04-28 00:20:00,0.492,0.208,236.240,2.200,0.000,0.000,0.0,8.2
2007-04-28 00:21:00,,,,,,,,
2007-04-28 00:22:00,,,,,,,,
2007-04-28 00:23:00,,,,,,,,
2007-04-28 00:24:00,,,,,,,,
2007-04-28 00:25:00,,,,,,,,
...

现在我们已经清理了数据集版本，我们可以使用可视化进一步调查它。

随着时间的推移观察模式

数据是多变量时间序列，理解时间序列的最佳方法是创建线图。

我们可以从为八个变量中的每一个创建单独的线图开始。

下面列出了完整的示例。

# line plots
from pandas import read_csv
from matplotlib import pyplot
# load the new file
dataset = read_csv('household_power_consumption.csv', header=0, infer_datetime_format=True, parse_dates=['datetime'], index_col=['datetime'])
# line plot for each variable
pyplot.figure()
for i in range(len(dataset.columns)):
	pyplot.subplot(len(dataset.columns), 1, i+1)
	name = dataset.columns[i]
	pyplot.plot(dataset[name])
	pyplot.title(name, y=0)
pyplot.show()

运行该示例将创建一个包含八个子图的单个图像，每个图对应一个变量。

这给了我们四分之一分钟观测的真正高水平。我们可以看到'Sub_metering_3'（环境控制）中可能没有直接映射到炎热或寒冷年份的有趣事情。也许安装了新系统。

有趣的是，'sub_metering_4'的贡献似乎随着时间的推移而减少，或呈现下降趋势，可能与'Sub_metering_3系列末尾的稳固增长相匹配”。

这些观察确实强调了在拟合和评估任何模型时遵守该数据的子序列的时间顺序的需要。

我们或许可以在'Global_active_power'和其他一些变量中看到季节性影响的波动。

有一些尖刻的用法可能与特定时期相匹配，例如周末。

功耗数据集中每个变量的线图

让我们放大并专注于'Global_active_power'或'_ 有功功率 _'。

我们可以为每年创建一个新的有效功率图，以查看这些年来是否存在任何共同模式。 2006 年的第一年，有不到一个月的数据，所以将其从情节中删除。

下面列出了完整的示例。

# yearly line plots
from pandas import read_csv
from matplotlib import pyplot
# load the new file
dataset = read_csv('household_power_consumption.csv', header=0, infer_datetime_format=True, parse_dates=['datetime'], index_col=['datetime'])
# plot active power for each year
years = ['2007', '2008', '2009', '2010']
pyplot.figure()
for i in range(len(years)):
	# prepare subplot
	ax = pyplot.subplot(len(years), 1, i+1)
	# determine the year to plot
	year = years[i]
	# get all observations for the year
	result = dataset[str(year)]
	# plot the active power for the year
	pyplot.plot(result['Global_active_power'])
	# add a title to the subplot
	pyplot.title(str(year), y=0, loc='left')
pyplot.show()

运行该示例将创建一个包含四个线图的单个图像，一个数据集中的每年全年（或大部分为全年）数据。

我们可以看到多年来的一些共同的总体模式，例如 2 月至 3 月左右和 8 月至 9 月左右，我们看到消费明显减少。

在夏季月份（北半球的年中），我们似乎也看到了下降的趋势，并且可能在冬季月份向地块的边缘消耗更多。这些可能显示出消费的年度季节性模式。

我们还可以在至少第一，第三和第四个图中看到一些缺失数据。

大多数年份的有功功率线图

我们可以继续放大消费量，并在 2007 年的 12 个月中查看有功功率。

这可能有助于梳理整个月的总体结构，例如每日和每周模式。

下面列出了完整的示例。

# monthly line plots
from pandas import read_csv
from matplotlib import pyplot
# load the new file
dataset = read_csv('household_power_consumption.csv', header=0, infer_datetime_format=True, parse_dates=['datetime'], index_col=['datetime'])
# plot active power for each year
months = [x for x in range(1, 13)]
pyplot.figure()
for i in range(len(months)):
	# prepare subplot
	ax = pyplot.subplot(len(months), 1, i+1)
	# determine the month to plot
	month = '2007-' + str(months[i])
	# get all observations for the month
	result = dataset[month]
	# plot the active power for the month
	pyplot.plot(result['Global_active_power'])
	# add a title to the subplot
	pyplot.title(month, y=0, loc='left')
pyplot.show()

运行该示例将创建一个包含 12 个线图的单个图像，2007 年每个月一个。

我们可以看到每个月内的日耗电的符号波。这很好，因为我们期望在功耗方面有某种日常模式。

我们可以看到，有很少的日子消费很少，例如 8 月和 4 月。这些可能代表住宅无人居住且耗电量最小的假期。

一年内所有月的有功功率线图

最后，我们可以放大一个级别，并仔细查看每日级别的功耗。

我们预计每天会有一些消费模式，也许一周内的天数差异。

下面列出了完整的示例。

# daily line plots
from pandas import read_csv
from matplotlib import pyplot
# load the new file
dataset = read_csv('household_power_consumption.csv', header=0, infer_datetime_format=True, parse_dates=['datetime'], index_col=['datetime'])
# plot active power for each year
days = [x for x in range(1, 20)]
pyplot.figure()
for i in range(len(days)):
	# prepare subplot
	ax = pyplot.subplot(len(days), 1, i+1)
	# determine the day to plot
	day = '2007-01-' + str(days[i])
	# get all observations for the day
	result = dataset[day]
	# plot the active power for the day
	pyplot.plot(result['Global_active_power'])
	# add a title to the subplot
	pyplot.title(day, y=0, loc='left')
pyplot.show()

运行该示例将创建一个包含 20 个线图的单个图像，一个用于 2007 年 1 月的前 20 天。

这些日子有共同之处;例如，很多天消费开始于凌晨 6 点到 7 点左右。

有些日子显示当天中午消费量下降，如果大多数人都不在家，这可能是有意义的。

我们确实看到有些日子有一些强烈的隔夜消费，在北半球，1 月可能与使用的供暖系统相匹配。

如预期的那样，一年中的时间，特别是它带来的季节和天气，将是对这些数据进行建模的重要因素。

一个月内 20 天的有功功率线图

时间序列数据分布

另一个需要考虑的重要领域是变量的分布。

例如，知道观测的分布是高斯分布还是其他分布可能是有趣的。

我们可以通过查看直方图来调查数据的分布。

我们可以通过为时间序列中的每个变量创建直方图来开始。

下面列出了完整的示例。

# histogram plots
from pandas import read_csv
from matplotlib import pyplot
# load the new file
dataset = read_csv('household_power_consumption.csv', header=0, infer_datetime_format=True, parse_dates=['datetime'], index_col=['datetime'])
# histogram plot for each variable
pyplot.figure()
for i in range(len(dataset.columns)):
	pyplot.subplot(len(dataset.columns), 1, i+1)
	name = dataset.columns[i]
	dataset[name].hist(bins=100)
	pyplot.title(name, y=0)
pyplot.show()

运行该示例会为 8 个变量中的每个变量创建一个单独的直方图。

我们可以看到，有功和无功功率，强度以及分计量功率都是偏向低瓦特小时或千瓦值的分布。

我们还可以看到电压数据的分布是强高斯分布的。

功耗数据集中每个变量的直方图

有功功率的分布似乎是双模态的，这意味着它看起来像有两组平均观察结果。

我们可以通过查看四年全年数据的有功功耗分布来进一步研究。

下面列出了完整的示例。

# yearly histogram plots
from pandas import read_csv
from matplotlib import pyplot
# load the new file
dataset = read_csv('household_power_consumption.csv', header=0, infer_datetime_format=True, parse_dates=['datetime'], index_col=['datetime'])
# plot active power for each year
years = ['2007', '2008', '2009', '2010']
pyplot.figure()
for i in range(len(years)):
	# prepare subplot
	ax = pyplot.subplot(len(years), 1, i+1)
	# determine the year to plot
	year = years[i]
	# get all observations for the year
	result = dataset[str(year)]
	# plot the active power for the year
	result['Global_active_power'].hist(bins=100)
	# zoom in on the distribution
	ax.set_xlim(0, 5)
	# add a title to the subplot
	pyplot.title(str(year), y=0, loc='right')
pyplot.show()

运行该示例将创建一个包含四个数字的单个图，每个图表在 2007 年到 2010 年之间的每一年。

我们可以看到这些年来有功功率消耗的分布看起来非常相似。分布确实是双峰的，一个峰值约为 0.3 千瓦，也许另一个峰值约为 1.3 千瓦。

分配到更高的千瓦值有一个长尾。它可能为将数据离散化并将其分为峰 1，峰 2 或长尾的概念敞开大门。这些用于一天或一小时的组或群集可能有助于开发预测模型。

大多数年份的有功功率直方图

所确定的群体可能在一年中的季节中变化。

我们可以通过查看一年中每个月的有功功率分布来研究这一点。

下面列出了完整的示例。

# monthly histogram plots
from pandas import read_csv
from matplotlib import pyplot
# load the new file
dataset = read_csv('household_power_consumption.csv', header=0, infer_datetime_format=True, parse_dates=['datetime'], index_col=['datetime'])
# plot active power for each year
months = [x for x in range(1, 13)]
pyplot.figure()
for i in range(len(months)):
	# prepare subplot
	ax = pyplot.subplot(len(months), 1, i+1)
	# determine the month to plot
	month = '2007-' + str(months[i])
	# get all observations for the month
	result = dataset[month]
	# plot the active power for the month
	result['Global_active_power'].hist(bins=100)
	# zoom in on the distribution
	ax.set_xlim(0, 5)
	# add a title to the subplot
	pyplot.title(month, y=0, loc='right')
pyplot.show()

运行该示例将创建一个包含 12 个图的图像，2007 年每个月一个。

我们每个月可以看到相同的数据分布。图中的轴似乎对齐（给定相似的比例），我们可以看到峰值在北半球温暖的月份向下移动，并在较冷的月份向上移动。

在 12 月至 3 月的较冷月份，我们还可以看到更大或更突出的尾部朝向更大的千瓦值。

一年内所有月份的有功功率直方图

关于建模的想法

现在我们知道了如何加载和探索数据集，我们可以提出一些关于如何建模数据集的想法。

在本节中，我们将在处理数据时仔细研究三个主要方面;他们是：

• 问题框架
• 数据准备
• 建模方法

问题框架

似乎没有关于数据集的开创性出版物来演示在预测性建模问题中构建数据的预期方法。

因此，我们可能会猜测可能使用这些数据的有用方法。

这些数据仅适用于单个家庭，但也许有效的建模方法可以推广到类似的家庭。

也许数据集最有用的框架是预测未来有效功耗的间隔。

四个例子包括：

• 预测第二天的每小时消耗量。
• 预测下周的每日消费量。
• 预测下个月的每日消费量。
• 预测下一年的月消费量。

通常，这些类型的预测问题称为多步预测。使用所有变量的模型可称为多变量多步预测模型。

这些模型中的每一个都不限于预测微小数据，而是可以将问题建模为所选预测分辨率或低于所选预测分辨率。

按规模预测消费可以帮助公用事业公司预测需求，这是一个广泛研究和重要的问题。

数据准备

在为建模准备这些数据时有很大的灵活性。

具体的数据准备方法及其益处实际上取决于所选择的问题框架和建模方法。不过，下面列出了可能有用的一般数据准备方法：

• 每日差异可用于调整数据中的每日循环。
• 年度差异可用于调整数据中的任何年度周期。
• 归一化可以帮助将具有不同单位的变量减少到相同的比例。

有许多简单的人为因素可能有助于数据的工程特征，反过来可能使特定的日子更容易预测。

一些例子包括：

• 指示一天中的时间，以说明人们回家的可能性。
• 指示一天是工作日还是周末。
• 指示某一天是否是北美公众假期。

这些因素对于预测月度数据可能要少得多，也许在每周数据的程度上要少得多。

更一般的功能可能包括：

• 指示季节，这可能导致使用的环境控制系统的类型或数量。

建模方法

对于这个问题，可能有四类方法可能很有趣;他们是：

• 朴素的方法。
• 经典线性方法。
• 机器学习方法。
• 深度学习方法。

朴素的方法

朴素的方法将包括做出非常简单但通常非常有效的假设的方法。

一些例子包括：

• 明天将和今天一样。
• 明天将与去年的这一天相同。
• 明天将是过去几天的平均值。

经典线性方法

经典线性方法包括对单变量时间序列预测非常有效的技术。

两个重要的例子包括：

• SARIMA
• ETS（三指数平滑）

他们需要丢弃其他变量，并将模型的参数配置或调整到数据集的特定框架。还可以直接支持与调整日常和季节性结构数据相关的问题。

机器学习方法

机器学习方法要求将问题构成监督学习问题。

这将要求将系列的滞后观察框架化为输入特征，从而丢弃数据中的时间关系。

可以探索一套非线性和集合方法，包括：

• k-最近邻居。
• 支持向量机
• 决策树
• 随机森林
• 梯度增压机

需要特别注意确保这些模型的拟合和评估保留了数据中的时间结构。这很重要，因此该方法无法通过利用未来的观测结果来“欺骗”。

这些方法通常与大量变量无关，可能有助于弄清楚是否可以利用其他变量并为预测模型增加价值。

深度学习方法

通常，神经网络在自回归类型问题上未被证明非常有效。

然而，诸如卷积神经网络的技术能够从原始数据（包括一维信号数据）自动学习复杂特征。并且诸如长短期记忆网络之类的循环神经网络能够直接学习输入数据的多个并行序列。

此外，这些方法的组合，例如 CNN LSTM 和 ConvLSTM，已经证明在时间序列分类任务上是有效的。

这些方法可能能够利用大量基于分钟的数据和多个输入变量。

进一步阅读

如果您希望深入了解，本节将提供有关该主题的更多资源。

• 个人家庭用电量数据集，UCI 机器学习库。
• 交流电源，维基百科。
• pandas.read_csv API

摘要

在本教程中，您发现了用于多步时间序列预测的家庭功耗数据集，以及如何使用探索性分析更好地理解原始数据。

具体来说，你学到了：

• 家庭用电量数据集，描述四年内单个房屋的用电量。
• 如何使用一系列线图来探索和理解数据集，用于数据分布的系列数据和直方图。
• 如何使用对问题的新理解来考虑预测问题的不同框架，可以准备数据的方式以及可以使用的建模方法。

你有任何问题吗？ 在下面的评论中提出您的问题，我会尽力回答。