时间序列预测——线性回归（上下界、异常检测），异常检测时候历史数据的输入选择是关键，使用过去历史值增加模型健壮性
data download:
https://github.com/nicolasmiller/pyculiarity/blob/master/tests/raw_data.csv
数据集样子：
y timestamp 1980-09-25 14:01:00 182.478 1980-09-25 14:02:00 176.231 1980-09-25 14:03:00 183.917 1980-09-25 14:04:00 177.798 1980-09-25 14:05:00 165.469 1980-09-25 14:06:00 181.878 1980-09-25 14:07:00 184.502 1980-09-25 14:08:00 183.303 1980-09-25 14:09:00 177.578 1980-09-25 14:10:00 171.641 1980-09-25 14:11:00 191.014 1980-09-25 14:12:00 184.068 1980-09-25 14:13:00 188.457 1980-09-25 14:14:00 175.739 1980-09-25 14:15:00 175.524 1980-09-25 14:16:00 189.128 1980-09-25 14:17:00 176.885 1980-09-25 14:18:00 167.140 1980-09-25 14:19:00 173.723 1980-09-25 14:20:00 168.460 1980-09-25 14:21:00 177.623 1980-09-25 14:22:00 183.888
做了shift处理前后：
y
timestamp
1980-09-25 14:01:00 182.478
1980-09-25 14:02:00 176.231
1980-09-25 14:03:00 183.917
1980-09-25 14:04:00 177.798
1980-09-25 14:05:00 165.469
1980-09-25 14:06:00 181.878
1980-09-25 14:07:00 184.502
1980-09-25 14:08:00 183.303
1980-09-25 14:09:00 177.578
1980-09-25 14:10:00 171.641
1980-09-25 14:11:00 191.014
1980-09-25 14:12:00 184.068
1980-09-25 14:13:00 188.457
1980-09-25 14:14:00 175.739
1980-09-25 14:15:00 175.524
1980-09-25 14:16:00 189.128
1980-09-25 14:17:00 176.885
1980-09-25 14:18:00 167.140
1980-09-25 14:19:00 173.723
1980-09-25 14:20:00 168.460
1980-09-25 14:21:00 177.623
1980-09-25 14:22:00 183.888
1980-09-25 14:23:00 167.487
1980-09-25 14:24:00 165.572
1980-09-25 14:25:00 170.480
1980-09-25 14:26:00 172.474
1980-09-25 14:27:00 166.448
1980-09-25 14:28:00 163.098
1980-09-25 14:29:00 163.544
1980-09-25 14:30:00 163.816
y lag_6 ... lag_23 lag_24
timestamp ...
1980-09-25 14:01:00 182.478 NaN ... NaN NaN
1980-09-25 14:02:00 176.231 NaN ... NaN NaN
1980-09-25 14:03:00 183.917 NaN ... NaN NaN
1980-09-25 14:04:00 177.798 NaN ... NaN NaN
1980-09-25 14:05:00 165.469 NaN ... NaN NaN
1980-09-25 14:06:00 181.878 NaN ... NaN NaN
1980-09-25 14:07:00 184.502 182.478 ... NaN NaN
1980-09-25 14:08:00 183.303 176.231 ... NaN NaN
1980-09-25 14:09:00 177.578 183.917 ... NaN NaN
1980-09-25 14:10:00 171.641 177.798 ... NaN NaN
1980-09-25 14:11:00 191.014 165.469 ... NaN NaN
1980-09-25 14:12:00 184.068 181.878 ... NaN NaN
1980-09-25 14:13:00 188.457 184.502 ... NaN NaN
1980-09-25 14:14:00 175.739 183.303 ... NaN NaN
1980-09-25 14:15:00 175.524 177.578 ... NaN NaN
1980-09-25 14:16:00 189.128 171.641 ... NaN NaN
1980-09-25 14:17:00 176.885 191.014 ... NaN NaN
1980-09-25 14:18:00 167.140 184.068 ... NaN NaN
1980-09-25 14:19:00 173.723 188.457 ... NaN NaN
1980-09-25 14:20:00 168.460 175.739 ... NaN NaN
1980-09-25 14:21:00 177.623 175.524 ... NaN NaN
1980-09-25 14:22:00 183.888 189.128 ... NaN NaN
1980-09-25 14:23:00 167.487 176.885 ... NaN NaN
1980-09-25 14:24:00 165.572 167.140 ... 182.478 NaN
1980-09-25 14:25:00 170.480 173.723 ... 176.231 182.478
1980-09-25 14:26:00 172.474 168.460 ... 183.917 176.231
1980-09-25 14:27:00 166.448 177.623 ... 177.798 183.917
1980-09-25 14:28:00 163.098 183.888 ... 165.469 177.798
1980-09-25 14:29:00 163.544 167.487 ... 181.878 165.469
1980-09-25 14:30:00 163.816 165.572 ... 184.502 181.878
代码：
# coding: utf-8
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from sklearn.model_selection import TimeSeriesSplit # you have everything done for you
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LassoCV, RidgeCV
# for time-series cross-validation set 5 folds
tscv = TimeSeriesSplit(n_splits=5)
def mean_absolute_percentage_error(y_true, y_pred):
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
def timeseries_train_test_split(X, y, test_size):
"""
Perform train-test split with respect to time series structure
"""
# get the index after which test set starts
test_index = int(len(X) * (1 - test_size))
X_train = X.iloc[:test_index]
y_train = y.iloc[:test_index]
X_test = X.iloc[test_index:]
y_test = y.iloc[test_index:]
return X_train, X_test, y_train, y_test
def plotModelResults(model, X_train, X_test, y_train, y_test, plot_intervals=False, plot_anomalies=False):
"""
Plots modelled vs fact values, prediction intervals and anomalies
"""
prediction = model.predict(X_test)
plt.figure(figsize=(15, 7))
plt.plot(prediction, "g", label="prediction", linewidth=2.0)
plt.plot(y_test.values, label="actual", linewidth=2.0)
if plot_intervals:
cv = cross_val_score(model, X_train, y_train,
cv=tscv,
scoring="neg_mean_absolute_error")
mae = cv.mean() * (-1)
deviation = cv.std()
scale = 20
lower = prediction - (mae + scale * deviation)
upper = prediction + (mae + scale * deviation)
plt.plot(lower, "r--", label="upper bond / lower bond", alpha=0.5)
plt.plot(upper, "r--", alpha=0.5)
if plot_anomalies:
anomalies = np.array([np.NaN] * len(y_test))
anomalies[y_test < lower] = y_test[y_test < lower]
anomalies[y_test > upper] = y_test[y_test > upper]
plt.plot(anomalies, "o", markersize=10, label="Anomalies")
error = mean_absolute_percentage_error(prediction, y_test)
plt.title("Mean absolute percentage error {0:.2f}%".format(error))
plt.legend(loc="best")
plt.tight_layout()
plt.grid(True);
plt.savefig("linear.png")
def plotCoefficients(model, X_train):
"""
Plots sorted coefficient values of the model
"""
coefs = pd.DataFrame(model.coef_, X_train.columns)
coefs.columns = ["coef"]
coefs["abs"] = coefs.coef.apply(np.abs)
coefs = coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)
plt.figure(figsize=(20, 12))
coefs.coef.plot(kind='bar')
plt.grid(True, axis='y')
plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles='dashed')
plt.savefig("linear-cov.png")
def code_mean(data, cat_feature, real_feature):
"""
Returns a dictionary where keys are unique categories of the cat_feature,
and values are means over real_feature
"""
return dict(data.groupby(cat_feature)[real_feature].mean())
def prepareData(series, lag_start, lag_end, test_size, target_encoding=False):
"""
series: pd.DataFrame
dataframe with timeseries
lag_start: int
initial step back in time to slice target variable
example - lag_start = 1 means that the model
will see yesterday's values to predict today
lag_end: int
final step back in time to slice target variable
example - lag_end = 4 means that the model
will see up to 4 days back in time to predict today
test_size: float
size of the test dataset after train/test split as percentage of dataset
target_encoding: boolean
if True - add target averages to the dataset
"""
# copy of the initial dataset
data = pd.DataFrame(series.copy())
data.columns = ["y"]
# lags of series
for i in range(lag_start, lag_end):
data["lag_{}".format(i)] = data.y.shift(i)
# datetime features
# data.index = data.index.to_datetime()
data["hour"] = data.index.hour
data["weekday"] = data.index.weekday
data['is_weekend'] = data.weekday.isin([5, 6]) * 1
if target_encoding:
# calculate averages on train set only
test_index = int(len(data.dropna()) * (1 - test_size))
data['weekday_average'] = list(map(
code_mean(data[:test_index], 'weekday', "y").get, data.weekday))
data["hour_average"] = list(map(
code_mean(data[:test_index], 'hour', "y").get, data.hour))
# drop encoded variables
# data.drop(["hour", "weekday"], axis=1, inplace=True)
# train-test split
y = data.dropna().y
X = data.dropna().drop(['y'], axis=1)
X_train, X_test, y_train, y_test = \
timeseries_train_test_split(X, y, test_size=test_size)
return X_train, X_test, y_train, y_test
def plt_linear():
data = pd.read_csv('raw_data.csv', usecols=['timestamp', 'count'])
data['timestamp'] = pd.to_datetime(data['timestamp'])
data.set_index("timestamp", drop=True, inplace=True)
data.rename(columns={'count': 'y'}, inplace=True)
X_train, X_test, y_train, y_test = \
prepareData(data, lag_start=6, lag_end=25, test_size=0.3, target_encoding=True)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# lr = LinearRegression()
lr = LassoCV(cv=tscv)
# lr = RidgeCV(cv=tscv)
# lr.fit(X_train_scaled, y_train)
"""
from xgboost import XGBRegressor
lr = XGBRegressor()
# lr = xgb.XGBRegressor(max_depth=5, learning_rate=0.1, n_estimators=160, silent=False, objective='reg:gamma')
"""
lr.fit(X_train_scaled, y_train)
"""
IMPORTANT
Generally tree-based models poorly handle trends in data, compared to linear models,
so you have to detrend your series first or use some tricks to make the magic happen.
Ideally make the series stationary and then use XGBoost, for example, you can forecast
trend separately with a linear model and then add predictions from xgboost to get final forecast.
"""
plotModelResults(lr, X_train=X_train_scaled, X_test=X_test_scaled, y_train=y_train, y_test=y_test, plot_intervals=True, plot_anomalies=True)
plotCoefficients(lr, X_train=X_train)
plt_linear()
可以看到相关系数！
重构代码，使其可以预测未来：
# coding: utf-8
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from sklearn.model_selection import TimeSeriesSplit
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LassoCV, RidgeCV
def mean_absolute_percentage_error(y_true, y_pred):
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
def timeseries_train_test_split(X, y, test_size, predict_size):
"""
Perform train-test split with respect to time series structure
"""
total_size = len(X)
# get the index after which test set starts
test_index = int(total_size * (1 - test_size))
X_train = X.iloc[:test_index]
y_train = y.iloc[:test_index]
X_test = X.iloc[test_index:total_size-predict_size]
y_test = y.iloc[test_index:total_size-predict_size]
X_predict = X.iloc[-predict_size:]
y_predict = y.iloc[-predict_size:]
return X_train, X_test, y_train, y_test, X_predict, y_predict
def predict_future(lr, X_predict, y_predict, lag_start, lag_end, scaler):
# for predict
# not OK for abnormal real value
y_predict[0:lag_start] = lr.predict(scaler.transform(X_predict.iloc[0:lag_start]))
# last_line = X_test.iloc[-1]
for i in range(lag_start, len(X_predict)):
# for i in range(0, len(X_predict)):
last_line = X_predict.iloc[i-1]
index = X_predict.index[i]
for j in range(lag_end-1, lag_start):
X_predict.at[index, "lag_{}".format(j)] = last_line["lag_{}".format(j-1)]
X_predict.at[index, "lag_{}".format(lag_start)] = y_predict[i-1]
y_predict[i] = lr.predict(scaler.transform([X_predict.iloc[i]]))[0]
return y_predict
def plot_results(y_predict, y, intervals, img_filename, plot_intervals=False, plot_anomalies=False, extra_plot=None):
"""
Plots modelled vs fact values, prediction intervals and anomalies
"""
assert len(y_predict) == len(y)
plt.figure(figsize=(15, 7))
# plt.plot(y.index, y_predict, "g", label="prediction", linewidth=3.0)
# plt.plot(y.index, y.values, label="actual", linewidth=1.0)
plt.plot(y.index, y_predict, ls='-', c='#0072B2', label='predicted y')
plt.plot(y.index, y.values, 'k.', label='y')
if extra_plot is not None:
# plt.plot(extra_plot.index, extra_plot.values, "y", label="future predict", linewidth=3.0)
plt.plot(extra_plot.index, extra_plot.values, 'y', label='predicted y')
if plot_intervals:
lower = y_predict - intervals
upper = y_predict + intervals
# plt.plot(y.index, lower, "r--", label="upper bond / lower bond", alpha=0.5)
# plt.plot(y.index, upper, "r--", alpha=0.5)
plt.fill_between(y.index, lower, upper, color='#0072B2', alpha=0.2, label='predicted upper/lower y')
if extra_plot is not None:
# plt.plot(extra_plot.index, extra_plot.values-intervals, "r--", label="upper bond / lower bond", alpha=0.5)
# plt.plot(extra_plot.index, extra_plot.values+intervals, "r--", alpha=0.5)
plt.fill_between(extra_plot.index, extra_plot.values-intervals, extra_plot.values+intervals,
color='#0072B2', alpha=0.2, label='predicted upper/lower y')
if plot_anomalies:
anomalies_lower = y[y < lower]
anomalies_upper = y[y > upper]
# plt.plot(anomalies_lower.index, anomalies_lower.values, "ro", markersize=10, label="Anomalies(+)")
# plt.plot(anomalies_upper.index, anomalies_upper.values, "ro", markersize=10, label="Anomalies(-)")
plt.plot(anomalies_lower.index, anomalies_lower.values, "rX", label='abnormal points')
plt.plot(anomalies_upper.index, anomalies_upper.values, "rX")
error = mean_absolute_percentage_error(y_predict, y)
plt.title("Mean absolute percentage error {0:.2f}%".format(error))
plt.legend(loc="best")
plt.tight_layout()
plt.grid(True);
plt.savefig(img_filename)
def plot_arg_importance(model, X_train):
"""
Plots sorted coefficient values of the model
"""
coefs = pd.DataFrame(model.coef_, X_train.columns)
coefs.columns = ["coef"]
coefs["abs"] = coefs.coef.apply(np.abs)
coefs = coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)
plt.figure(figsize=(20, 12))
coefs.coef.plot(kind='bar')
plt.grid(True, axis='y')
plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles='dashed')
plt.savefig("linear-cov.png")
def code_mean(data, cat_feature, real_feature):
"""
Returns a dictionary where keys are unique categories of the cat_feature,
and values are means over real_feature
"""
return dict(data.groupby(cat_feature)[real_feature].mean())
def prepare_data(series, lag_start, lag_end, test_size, target_encoding=False, days_to_predict=1):
"""
series: pd.DataFrame
dataframe with timeseries
lag_start: int
initial step back in time to slice target variable
example - lag_start = 1 means that the model
will see yesterday's values to predict today
lag_end: int
final step back in time to slice target variable
example - lag_end = 4 means that the model
will see up to 4 days back in time to predict today
test_size: float
size of the test dataset after train/test split as percentage of dataset
target_encoding: boolean
if True - add target averages to the dataset
"""
last_date = series["timestamp"].max()
def make_future_date(periods, freq='D'):
"""Simulate the trend using the extrapolated generative model.
Parameters
----------
periods: Int number of periods to forecast forward.
freq: Any valid frequency for pd.date_range, such as 'D' or 'M'.
Returns
-------
pd.Dataframe that extends forward from the end of self.history for the
requested number of periods.
"""
dates = pd.date_range(
start=last_date,
periods=periods + 1, # An extra in case we include start
freq=freq)
dates = dates[dates > last_date] # Drop start if equals last_date
return dates[:periods] # Return correct number of periods
predict_points = days_to_predict * 1440 # 1 day = 60*24 minutes
future_dates = make_future_date(periods=predict_points, freq='T')
df_future = pd.DataFrame({"timestamp": future_dates, "y": np.zeros(len(future_dates))})
data = pd.concat([series, df_future])
data.set_index("timestamp", drop=True, inplace=True)
# data = pd.DataFrame(series.copy())
# data.columns = ["y"]
print(data[:30])
# lags of series
for i in range(lag_start, lag_end):
data["lag_{}".format(i)] = data.y.shift(i)
print(data[:30])
# datetime features
# data.index = data.index.to_datetime()
data["hour"] = data.index.hour
data["weekday"] = data.index.weekday
data['is_weekend'] = data.weekday.isin([5, 6]) * 1
if target_encoding:
# calculate averages on train set only
test_index = int(len(data.dropna()) * (1 - test_size))
data['weekday_average'] = list(map(
code_mean(data[:test_index], 'weekday', "y").get, data.weekday))
data["hour_average"] = list(map(
code_mean(data[:test_index], 'hour', "y").get, data.hour))
# drop encoded variables
# data.drop(["hour", "weekday"], axis=1, inplace=True)
# train-test split
y = data.dropna().y
X = data.dropna().drop(['y'], axis=1)
X_train, X_test, y_train, y_test, X_predict, y_predict = \
timeseries_train_test_split(X, y, test_size=test_size, predict_size=predict_points)
return X_train, X_test, y_train, y_test, X_predict, y_predict
def calculate_intevals(lr, X_train, y_train, tscv):
cv = cross_val_score(lr, X_train, y_train,
cv=tscv,
scoring="neg_mean_absolute_error")
mae = cv.mean() * (-1)
deviation = cv.std()
scale = 30
return mae + scale * deviation
def plt_linear():
data = pd.read_csv('raw_data.csv', usecols=['timestamp', 'count'])
# input format
data['timestamp'] = pd.to_datetime(data['timestamp'])
data = data.sort_values('timestamp')
data.rename(columns={'count': 'y'}, inplace=True)
lag_start = 1
lag_end = 100
X_train, X_test, y_train, y_test, X_predict, y_predict = \
prepare_data(data, lag_start=lag_start, lag_end=lag_end, test_size=0.3, target_encoding=True)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# for time-series cross-validation set 5 folds
tscv = TimeSeriesSplit(n_splits=5)
# lr = LinearRegression()
lr = LassoCV(cv=tscv)
# lr = RidgeCV(cv=tscv)
lr.fit(X_train_scaled, y_train)
intervals = calculate_intevals(lr, X_train, y_train, tscv)
"""
y_test_predict = lr.predict(X_test_scaled)
plot_results(y_predict=y_test_predict, y=y_test, intervals=intervals, img_filename="linear-test-result.png", plot_intervals=True, plot_anomalies=True)
"""
y2 = lr.predict(np.concatenate((X_train_scaled, X_test_scaled)))
y = pd.concat([y_train, y_test])
# plot_results(y_predict=y2, y=y, intervals=intervals, img_filename="linear-all-result.png", plot_intervals=True, plot_anomalies=True)
y_future = predict_future(lr, X_predict, y_predict, lag_start, lag_end, scaler)
plot_results(y_predict=y2, y=y, intervals=intervals, img_filename="linear-all.png", plot_intervals=True, plot_anomalies=True, extra_plot=y_future)
plot_arg_importance(lr, X_train=X_train)
plt_linear()
绘图：
尤其关键的是，lag_start, lag_end参数，如果lag_start=1，则表示使用前一个时刻输入，会导致模型过拟合，因此上面设置了lag_start=60表示使用一个小时前的数据来预测，防止过拟合。
来lag_start=1的效果：
可以看到，前一个时刻数据值的重要性！因此，最后做趋势预测的时候出现了重大失误，模型过拟合了。
bug修复：
# coding: utf-8
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from sklearn.model_selection import TimeSeriesSplit
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LassoCV, RidgeCV
from sklearn.ensemble import GradientBoostingRegressor
def mean_absolute_percentage_error(y_true, y_pred):
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
def predict_future(lr, X_predict, y_predict, lag_start, lag_end, scaler):
# for predict
y_predict[0:lag_start] = lr.predict(scaler.transform(X_predict.iloc[0:lag_start]))
for i in range(lag_start, len(X_predict)):
last_line = X_predict.iloc[i-1]
index = X_predict.index[i]
for j in range(lag_end-1, lag_start, -1):
X_predict.at[index, "lag_{}".format(j)] = last_line["lag_{}".format(j-1)]
X_predict.at[index, "lag_{}".format(lag_start)] = y_predict[i-1]
y_predict[i] = lr.predict(scaler.transform([X_predict.iloc[i]]))[0]
return y_predict
def plot_results(y_predict, y, intervals, img_filename, plot_intervals=False, plot_anomalies=False, extra_plot=None):
"""
Plots modelled vs fact values, prediction intervals and anomalies
"""
assert len(y_predict) == len(y)
plt.figure(figsize=(15, 7))
# plt.plot(y.index, y_predict, "g", label="prediction", linewidth=3.0)
# plt.plot(y.index, y.values, label="actual", linewidth=1.0)
plt.plot(y.index, y_predict, ls='-', c='#0072B2', label='predicted y')
plt.plot(y.index, y.values, 'k.', label='y')
if extra_plot is not None:
# plt.plot(extra_plot.index, extra_plot.values, "y", label="future predict", linewidth=3.0)
plt.plot(extra_plot.index, extra_plot.values, 'y', label='predicted y')
if plot_intervals:
lower = y_predict - intervals
upper = y_predict + intervals
# plt.plot(y.index, lower, "r--", label="upper bond / lower bond", alpha=0.5)
# plt.plot(y.index, upper, "r--", alpha=0.5)
plt.fill_between(y.index, lower, upper, color='#0072B2', alpha=0.2, label='predicted upper/lower y')
if extra_plot is not None:
# plt.plot(extra_plot.index, extra_plot.values-intervals, "r--", label="upper bond / lower bond", alpha=0.5)
# plt.plot(extra_plot.index, extra_plot.values+intervals, "r--", alpha=0.5)
plt.fill_between(extra_plot.index, extra_plot.values-intervals, extra_plot.values+intervals,
color='#0072B2', alpha=0.2, label='predicted upper/lower y')
if plot_anomalies:
anomalies_lower = y[y < lower]
anomalies_upper = y[y > upper]
# plt.plot(anomalies_lower.index, anomalies_lower.values, "ro", markersize=10, label="Anomalies(+)")
# plt.plot(anomalies_upper.index, anomalies_upper.values, "ro", markersize=10, label="Anomalies(-)")
plt.plot(anomalies_lower.index, anomalies_lower.values, "rX", label='abnormal points')
plt.plot(anomalies_upper.index, anomalies_upper.values, "rX")
error = mean_absolute_percentage_error(y_predict, y)
plt.title("Mean absolute percentage error {0:.2f}%".format(error))
plt.legend(loc="best")
plt.tight_layout()
plt.grid(True)
plt.savefig(img_filename)
def plot_arg_importance(model, X_train, img_filename="linear-cov.png"):
"""
Plots sorted coefficient values of the model
"""
coefs = pd.DataFrame(model.coef_, X_train.columns)
coefs.columns = ["coef"]
coefs["abs"] = coefs.coef.apply(np.abs)
coefs = coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)
plt.figure(figsize=(20, 12))
coefs.coef.plot(kind='bar')
plt.grid(True, axis='y')
plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles='dashed')
plt.savefig(img_filename)
def code_mean(data, cat_feature, real_feature):
"""
Returns a dictionary where keys are unique categories of the cat_feature,
and values are means over real_feature
"""
return dict(data.groupby(cat_feature)[real_feature].mean())
def prepare_data(series, lag_start, lag_end, test_size, target_encoding=False, days_to_predict=2):
"""
series: pd.DataFrame
dataframe with timeseries
lag_start: int
initial step back in time to slice target variable
example - lag_start = 1 means that the model
will see yesterday's values to predict today
lag_end: int
final step back in time to slice target variable
example - lag_end = 4 means that the model
will see up to 4 days back in time to predict today
test_size: float
size of the test dataset after train/test split as percentage of dataset
target_encoding: boolean
if True - add target averages to the dataset
"""
last_date = series["timestamp"].max()
def make_future_date(periods, freq='D'):
"""Simulate the trend using the extrapolated generative model.
Parameters
----------
periods: Int number of periods to forecast forward.
freq: Any valid frequency for pd.date_range, such as 'D' or 'M'.
Returns
-------
pd.Dataframe that extends forward from the end of self.history for the
requested number of periods.
"""
dates = pd.date_range(
start=last_date,
periods=periods + 1, # An extra in case we include start
freq=freq)
dates = dates[dates > last_date] # Drop start if equals last_date
return dates[:periods] # Return correct number of periods
predict_points = days_to_predict * 1440 # 1 day = 60*24 minutes
future_dates = make_future_date(periods=predict_points, freq='T')
df_future = pd.DataFrame({"timestamp": future_dates, "y": np.zeros(len(future_dates))})
data = pd.concat([series, df_future])
data.set_index("timestamp", drop=True, inplace=True)
# data = pd.DataFrame(series.copy())
# data.columns = ["y"]
# print(data[:30])
# lags of series
for i in range(lag_start, lag_end):
data["lag_{}".format(i)] = data.y.shift(i)
# print(data[:30])
# datetime features
# data.index = data.index.to_datetime()
data["hour"] = data.index.hour
data["weekday"] = data.index.weekday
data['is_weekend'] = data.weekday.isin([5, 6]) * 1
test_index = int(len(series) * (1 - test_size))
if target_encoding:
# calculate averages on train set only
data['weekday_average'] = list(map(
code_mean(data[:test_index], 'weekday', "y").get, data.weekday))
data["hour_average"] = list(map(
code_mean(data[:test_index], 'hour', "y").get, data.hour))
# drop encoded variables
# data.drop(["hour", "weekday"], axis=1, inplace=True)
# train-test split
y = data.dropna().y
X = data.dropna().drop(['y'], axis=1)
total_size = len(X)
# get the index after which test set starts
X_train = X.iloc[:test_index]
y_train = y.iloc[:test_index]
X_test = X.iloc[test_index:total_size-predict_points]
y_test = y.iloc[test_index:total_size-predict_points]
X_predict = X.iloc[-predict_points:]
y_predict = y.iloc[-predict_points:]
return X_train, X_test, y_train, y_test, X_predict, y_predict
def mean_absolute_error(y_true, y_pred):
abs_err = np.abs(y_true-y_pred)
return np.mean(abs_err), np.std(abs_err)
def calculate_intervals2(y_true, y_pred, scale):
mae, std = mean_absolute_error(y_true, y_pred)
return mae + scale * std
def calculate_intervals(lr, X_train, y_train, tscv, scale):
cv = cross_val_score(lr, X_train, y_train,
cv=tscv,
scoring="neg_mean_squared_error")
mae = cv.mean() * (-1)
deviation = cv.std()
return mae + scale * deviation
def linear_predict(data_frame, interval_scale=30, lag_start=60, lag_end=100, days_to_predict=2):
"""
predict time series data using linear model.
:param data_frame: input data frame. Must have timestamp and y columns. Also set index with timestamp.
:param interval_scale: interval range scale.
:param lag_start:
initial step back in time to slice target variable
example - lag_start = 1 means that the model
will see yesterday's values to predict today
:param lag_end:
final step back in time to slice target variable
example - lag_end = 4 means that the model
will see up to 4 days back in time to predict today
:param days_to_predict: predicted days in future/
:return:
history_predict: history predicted value (pandas.Series)
y_future: predicted value for future (pandas.Series)
intervals: upper and lower range (float).
"""
assert "timestamp" in data_frame.columns
assert "y" in data_frame.columns
X_train, X_test, y_train, y_test, X_predict, y_predict = \
prepare_data(data_frame, lag_start=lag_start, lag_end=lag_end, test_size=0.3, target_encoding=True, days_to_predict=days_to_predict)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# for time-series cross-validation set 5 folds
tscv = TimeSeriesSplit(n_splits=5)
# lr = LinearRegression()
lr = GradientBoostingRegressor(n_estimators=100, learning_rate=0.1, max_depth=1, random_state=666)
# lr = LassoCV(cv=tscv)
# lr = RidgeCV(cv=tscv)
lr.fit(X_train_scaled, y_train)
# plot_arg_importance(lr, X_train=X_train, img_filename="linear-cov.png")
y_future = predict_future(lr, X_predict, y_predict, lag_start, lag_end, scaler)
y = lr.predict(np.concatenate((X_train_scaled, X_test_scaled)))
y_history = pd.concat([y_train, y_test])
# intervals = calculate_intervals(lr, X_train, y_train, tscv, scale=interval_scale)
intervals = calculate_intervals2(y_history, y, interval_scale)
# anomalies_lower = y_history[y_history<y-intervals]
# anomalies_upper = y_history[y_history>y+intervals]
assert len(y_history.index) == len(y)
return pd.Series(data=y, index=y_history.index, name="history_predict"), y_future, intervals
def get_anoms(y_real, y_predict, intervals, lower_ratio=1.0, upper_ratio=1.0):
"""
calculat anomal point using predicted value and interval
:param y_real: real value (pandas.Series)
:param y_predict: predicted value (pandas.Series)
:param intervals: upper and lower bound range (float)
:param lower_ratio: lower ratio you want to scale
:param upper_ratio: upper ratio you want scale
:return:
anoms_lower, anoms_upper (pandas.Series)
"""
anomalies_lower_index,anomalies_lower_val = [], []
anomalies_upper_index,anomalies_upper_val = [], []
for timestamp, expect_val in zip(y_predict.index, y_predict.values):
real_val = y_real.loc[timestamp]
if (expect_val-intervals) * lower_ratio > real_val:
anomalies_lower_index.append(timestamp)
anomalies_lower_val.append(real_val)
if (expect_val+intervals) * upper_ratio < real_val:
anomalies_upper_index.append(timestamp)
anomalies_upper_val.append(real_val)
return pd.Series(data=anomalies_lower_val, index=anomalies_lower_index),\
pd.Series(data=anomalies_upper_val, index=anomalies_upper_index)
def plot_history_and_future(y_predict, y_real, intervals, anomalies_lower, anomalies_upper, predicted_future, img_filename, need_lower=True):
"""
Plots modelled vs fact values, prediction intervals and anomalies
"""
assert len(y_predict) < len(y_real)
plt.figure(figsize=(15, 7))
plt.plot(y_predict.index, y_predict.values, ls='-', c='#0072B2', label='predicted y')
plt.plot(y_real.index, y_real.values, 'k.', label='y')
if need_lower:
plt.fill_between(y_predict.index, y_predict.values - intervals, y_predict.values + intervals, color='#0072B2', alpha=0.2, label='predicted upper/lower y')
else:
plt.fill_between(y_predict.index, 0, y_predict.values + intervals, color='#0072B2', alpha=0.2, label='predicted upper/lower y')
plt.plot(predicted_future.index, predicted_future.values, 'y', label='predicted y')
if need_lower:
plt.fill_between(predicted_future.index, predicted_future.values-intervals, predicted_future.values+intervals,
color='#0072B2', alpha=0.2)
else:
plt.fill_between(predicted_future.index, 0, predicted_future.values+intervals,
color='#0072B2', alpha=0.2)
if need_lower:
plt.plot(anomalies_lower.index, anomalies_lower.values, "rX", label='abnormal points')
plt.plot(anomalies_upper.index, anomalies_upper.values, "rX")
else:
plt.plot(anomalies_upper.index, anomalies_upper.values, "rX", label="abnormal points")
error = mean_absolute_percentage_error(y_predict, y_real)
plt.title("Mean absolute percentage error {0:.2f}%".format(error))
plt.legend(loc="best")
plt.tight_layout()
plt.grid(True)
plt.savefig(img_filename)
if __name__ == "__main__":
data = pd.read_csv('raw_data.csv', usecols=['timestamp', 'count'])
# input format
data['timestamp'] = pd.to_datetime(data['timestamp'])
data = data.sort_values('timestamp')
data.rename(columns={'count': 'y'}, inplace=True)
data.set_index("timestamp", drop=False, inplace=True)
y_predict, y_future, intervals = linear_predict(data, interval_scale=5, lag_start=60, lag_end=100, days_to_predict=3)
anomalies_lower, anomalies_upper = get_anoms(data['y'], y_predict, intervals)
plot_history_and_future(y_predict=y_predict, y_real=data['y'], intervals=intervals, anomalies_lower=anomalies_lower,
anomalies_upper=anomalies_upper, predicted_future=y_future, img_filename="linear.png", need_lower=True)
修复了-1的问题：
for j in range(lag_end-1, lag_start, -1):
lag的bug。
此外使用梯度提升树模型做回归，目前看效果略好于其他模型，线性回归模型很健壮，但是在特殊情况下会出现网络流预测值为为负数的情形，根因还没有找到，而梯度提升树没有这个问题，但是GBT在数据平稳，预测应该是常数的时候会出现上升情形（线性回归也有这个问题，WHY？？？）。如下图所示数据情形:
