大家好,我是yudengwu
时间序列特征构造
时间序列问题,首先不管是回归问题,还是分类问题。
一个模型的好坏,决定因素由数据集的大小,特征值的选取和处理,算法。
其中最重要的是特征值的选取和处理。
今天余总来讲解下时间序列的特征构造问题。
该特征构造部分可以用于其他数值数据。
时间序列特征构造分类为 :时间特征,时间历史特征,时间交叉特征
时间特征
连续时间:
持续时间,间隔时间
离散时间:
年,季度,季节,月,星期,日,等
节假日,节假日第几天
上午,早上,中午,晚上
年初,年末,月初,月末,周内,周末
是否高峰时段,是否上班
时间历史特征
统计值:
四分位数,中位数,平均值,偏度,峰度,离散系数
同期值:
如今天的8点和昨天的8点,
这周一和上周一
时间交叉特征
类别特征与类别特征:笛卡尔积
连续特征与类别特征:离散笛卡尔积,聚合特征(聚合指聚类后的类别特征)
连续特征与连续特征:一阶差分,二阶差分(即把数据和前后做差值,差值作为特征,一次差值为一阶差分)
时间特征构造代码
使用的包 pandas
个人pandas教程总结
pandas常用函数,个人常用的
时间特征构造
原始数据total_df.csv
给出的数据只有时间和负荷,
我们要从时间里提取出年,月,日,节假日,工作日等时间特征。
首先把文档里的时间数据格式由字符串转换到时间
total_df.loc[:, 'record_date'] = pd.to_datetime(total_df['record_date'],format='%Y-%m-%dT%H:%M:%S')
利用python里的time 函数模块,构造强时间指代特征:星期几(dow)、几号(dom)、几月(month)、哪一年(year)。代码如下:
# 几时
# 星期几,其中dow为构造的新列
total_df.loc[:, 'dow'] = total_df['record_date'].apply(lambda x: x.dayofweek)
# 几号,dom为构造的新列
total_df.loc[:, 'dom'] = total_df['record_date'].apply(lambda x: x.day)
# 几月,month为构造的新列
total_df.loc[:, 'month'] = total_df['record_date'].apply(lambda x: x.month)
# 几年,year为构造的新列
total_df.loc[:, 'year'] = total_df['record_date'].apply(lambda x: x.year)
利用time函数模块添加周末特征:通过0,1化添加周末特征,周末的特征为1,将文字数值化传入模型,具体代码如下:
#weekend,weekend_sat,weekend_sun 为构造的新列
total_df.loc[:,'weekend'] = 0
total_df.loc[:,'weekend_sat'] = 0
total_df.loc[:,'weekend_sun'] = 0
total_df.loc[(total_df['dow']>4), 'weekend'] = 1
total_df.loc[(total_df['dow']==5), 'weekend_sat'] = 1
total_df.loc[(total_df['dow']==6), 'weekend_sun'] = 1
添加特征:添加上下中旬,其特征分别表示为1/2/3,代码如下:
def period_of_month(day):
if day in range(1, 11):
return 1
if day in range(11, 21):
return 2
else:
return 3
#period_of_month为构造的新列
total_df.loc[:, 'period_of_month'] = total_df['dom'].apply(lambda x: period_of_month(x))
添加上半月、下半月时间特征,其特征分别表示为1/2,代码如下:
def period2_of_month(day):
if day in range(1, 16):
return 1
else:
return 2
#period2_of_month为构造的新列
total_df.loc[:, 'period2_of_month'] = total_df['dom'].apply(lambda x: period2_of_month(x))
利用python 里的time函数模块再自定义函数,添加国庆节假日特征,如果满足月份为10,日期小于8,则为国庆节。代码如下:
total_df.loc[:, 'festival'] = 0
total_df.loc[(total_df.month == 10) & (total_df.dom < 8), 'festival'] = 1
天气特征构造
天气啥的一般是文字,(如晴,多云…)文字程序无法直接输入,需要先One-hot化。
one-hot方法有很多,介绍下pandas 包的吧
介绍pd.get_dummies
示例
import pandas as pd
df = pd.DataFrame([
['green' , 'A'],
['red' , 'B'],
['blue' , 'A']])
df.columns = ['color', 'class']
print('处理前')
print(df)
df1=pd.get_dummies(df)
print('处理后')
print(df1)
一阶差分,二阶差分啥的,pandas shift,diff就能搞定
data['use1']=data['use'].shift(1)#向下平移1
data['use2']=data['use'].shift(2)#向下平移2
data['use3']=data['use'].shift(3)##向下平移3
data['use1差值']=data['use'].diff(1)#前后 t和t-1差值
data['use2差值']=data['use'].diff(2)##t和t-2差值
data= data.dropna()#去空值
统计特征啥的,pandas 有函数了调,不讲解,如果这都不会,自己拿豆腐撞墙吧。
聚合特征:首先把各个样本聚类,然后把类别信息作为特征之一,传入到输入中,作为输入的一部分。
例子:负荷预测
处理好后的数据样子截图
数据已经对天气,节假日进行了处理。
误差函数
#statistic.py
import numpy as np
import math
import scipy.stats as stats
#计算预测值真实值的均方误差
#actual 实际值
def meanSquareError(actual,pred):
if (not len(actual) == len(pred) or len(actual) == 0):
return -1.0
total = 0.0
for x in range(len(actual)):
total += math.pow(actual[x]-pred[x],2)
return total/len(actual)
# actual values实际值
#计算mse误差
def mse(actual,pred):
if (not len(actual) == len(pred) or len(actual) == 0):
return -1.0
total = 0.0
for x in range(len(actual)):
total += math.pow(actual[x]-pred[x],2)
return total/(len(actual)*1000000)
# 计算预测值与实际值之间的归一化均方根误差(NRMSE)
#
def normRmse(actual,pred):
if (not len(actual) == len(pred) or len(actual) == 0):
return -1.0
sumSquares = 0.0
maxY = actual[0]
minY = actual[0]
for x in range(len(actual)):
sumSquares += math.pow(pred[x]-actual[x],2.0)
maxY = max(maxY,actual[x])
minY = min(minY,actual[x])
return math.sqrt(sumSquares/len(actual))/(maxY-minY)
# 根据实际值计算预测值的均方根误差(RMSE)
def Rmse(actual,pred):
if (not len(actual) == len(pred) or len(actual) == 0):
return -1.0
sumSquares = 0.0
for x in range(len(actual)):
sumSquares += math.pow(pred[x]-actual[x],2.0)
return math.sqrt(sumSquares/len(actual))
#计算相对于实际值的预测值的绝对百分误差(MAPE)
def mape(actual,pred):
if (not len(actual) == len(pred) or len(actual) == 0):
return -1.0
total = 0.0
for x in range(len(actual)):
total += abs((actual[x]-pred[x])/actual[x])
return total/len(actual)
#计算相对于实际值的预测值的绝对百分误差(MAPE)
def mae(actual,pred):
if (not len(actual) == len(pred) or len(actual) == 0):
return -1.0
total = 0.0
for x in range(len(actual)):
total += abs(actual[x]-pred[x])
return total/len(actual)
预测代码
import math
from tools import statistics#引入误差函数
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
# 将时间序列转换为监督学习问题
#n_in=24时即为将24小时的所有数据去预测24小时后的数据。
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
n_vars = 1 if type(data) is list else data.shape[1]
df = DataFrame(data)
cols, names = list(), list()
# input sequence (t-n, ... t-1)
for i in range(n_in, 0, -1):
cols.append(df.shift(i))#当
names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]
# forecast sequence (t, t+1, ... t+n)
for i in range(0, n_out):
cols.append(df.shift(-i))#向上平移即将下一天当前时间的actual功率平移当当前天,作为输出标签
if i == 0:
names += [('var%d(t)' % (j+1)) for j in range(n_vars)]
else:
names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)]
# put it all together
agg = concat(cols, axis=1)
agg.columns = names
# drop rows with NaN values
if dropnan:
agg.dropna(inplace=True)
return agg
# 加载数据
dataset = read_csv('load.csv', header=0, index_col=0)
values = dataset.values
# 整数编码
# encoder = LabelEncoder()
# values[:,4] = encoder.fit_transform(values[:,4])
# ensure all data is float
values = values.astype('float32')
# normalize features 归一化特征
# scaler = MinMaxScaler(feature_range=(0, 1))
# scaled = scaler.fit_transform(values)
matrix_load = np.array(values)
matrix_load[:,0] = matrix_load[:,0]/1000#daily load 加载日负载数据,并将第一列功率每一个数除以1000
# 异常值处理
k = 0
for j in range(0, matrix_load.shape[0]):
##如何一个数与前一个相差2(即2000)和与后一个相差2(2000)
if(abs(matrix_load[j,0]-matrix_load[j-1,0])>2 and abs(matrix_load[j,0]-matrix_load[j+1,0])>2):
k = k + 1
##则这个数=(前一个+后一个)/2+(前一天的前一个+前一天的后一个)/2
matrix_load[j,0] = (matrix_load[j - 1,0] + matrix_load[j + 1,0]) / 2 + matrix_load[j - 24,0] - matrix_load[j - 24 - 1,0] / 2
sum = 0
num = 0
for t in range(1,8):
if(j - 24*t >= 0):
num = num + 1
sum = sum + matrix_load[j - 24*t,0]#前一周以内每日同时数据之和(如星期一的4点,星期二的4点...)
if((j + 24*t) < matrix_load.shape[0]):
num = num + 1
sum = sum + matrix_load[j + 24*t,0]#前一周的+后一周每日同时数据之和
sum = sum / num#前后一周每日同时数据平均值
# #如果当前数据与平均值相差超过3,则当前数据=平均值加或减3
if(abs(matrix_load[j,0] - sum)>3):
k = k + 1
if(matrix_load[j,0] > sum): matrix_load[j,0] = sum + 3
else: matrix_load[j,0] = sum - 3
print(k)
# print(matrix_load[1,0])
# shift all data by mean 去均值
shifted_value = matrix_load[:,0].mean()
matrix_load[:,0] -= shifted_value
for i in range(1, 9):
matrix_load[:, i] = matrix_load[:, i] / 10
shifted_valuei = matrix_load[:,i].mean()
matrix_load[:,i] -= shifted_valuei
print(matrix_load.shape)
# frame as supervised learning
before = 24
end = 1
reframed = series_to_supervised(matrix_load, before, end)
print(reframed.shape)
# drop columns we don't want to predict
for i in range(76):
reframed.drop(reframed.columns[[-1]], axis=1, inplace=True)
print(reframed.shape)
print(reframed.head())
# 划分训练集测试集
values = reframed.values
n_train_hours = values.shape[0] - 166*24
train = values[:n_train_hours, :]
test = values[n_train_hours:, :]
print('train:',n_train_hours,'test:',166*24)
# 划分输入和输出
train_X, train_y = train[:, :-end], train[:, -end]
print(train_X.shape, train_y.shape)
test_X, test_y = test[:, :-end], test[:, -end]
#改变数据形状 3D [samples, timesteps, features]
train_X = train_X.reshape((train_X.shape[0], before*77))
test_X = test_X.reshape((test_X.shape[0], before*77))
print(train_X.shape, train_y.shape, test_X.shape, test_y.shape)
# svr
# kernel='linear'时,为线性核,C越大分类效果越好,但有可能会过拟合(defaul C=1)
# kernel='rbf'时(default),为高斯核radial basis function,gamma值越小,分类界面越连续;gamma值越大,分类界面越“散”,分类效果越好,但有可能会过拟合
kernelList = ["rbf"]
names = ["true","radial basis"]
preds = []
preds.append(test_y)
for i in range(len(kernelList)):
clf = svm.SVR(C=2.0, kernel=kernelList[i])
clf.fit(train_X, train_y)
predicted_values = clf.predict(test_X)
mape = statistics.mape((test_y + shifted_value) * 1000, (predicted_values + shifted_value) * 1000)
print('MAPE is ', mape)
mae = statistics.mae((test_y + shifted_value) * 1000, (predicted_values + shifted_value) * 1000)
print('MAE is ', mae)
mse = statistics.meanSquareError((test_y + shifted_value) * 1000, (predicted_values + shifted_value) * 1000)
print('MSE is ', mse)
rmse = math.sqrt(mse)
print('RMSE is ', rmse)
nrmse = statistics.normRmse((test_y + shifted_value) * 1000, (predicted_values + shifted_value) * 1000)
print('NRMSE is ', nrmse)
preds.append(predicted_values)
# 结果
fig = plt.figure()
colors = ["g","r","b","c","m","y","k","w"]
legendVars = []
for j in range(len(preds)):
print(j)
x, = plt.plot(preds[j]+shifted_value, color=colors[j])
legendVars.append(x)
plt.title("Svm (rbf) Multi")
plt.xlabel('Time(hour)')
plt.ylabel('Electricity load (*1e3 kW)')
plt.legend(legendVars, names)
plt.show()
fig.savefig('svr_mul_result.jpg', bbox_inches='tight')
说明:
本负荷预测:在特征方面有
1.天气,节假日,周末,工作日
2.将前面24小时的负荷日期气温等因素 来预测第25小时的负荷。有一定的同期日特征,可以考虑其他同期日,如上一周的同日同时等,这个看自己取舍。
当然你也可以考虑其他特征,特征并不是越多越好。
时间特征不仅仅可以用在回归上,分类啥的也可以。
电气专业的计算机萌新,写博文不容易。如果你觉得本文对你有用,请点个赞支持下,谢谢。