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？？？）。如下圖所示數據情形:
