TensorFlow Probability (TFP)がリリースされてからしばらく経ちますが、最近になってこんなモジュールが公開されたと知りました。
Framework for Bayesian structural time series modelsと題されている通りで、ズバリTFPでベイズ構造時系列モデルを推定するためのモジュールです。
ベイズ構造時系列モデルの推定は以前このブログでもRStanやRの{bsts}パッケージを使って散々やりましたが、そもそもTFPはサンプリング手法にHMCを使うなど割とStanとの共通点が多く、またこのtfp.stsモジュール自体がそのクラスのネーミングからしてどう見ても{bsts}を意識したものなので、それらを時系列データ分析に使ってきた人にとっては取っ付きやすいかもしれません。
ということで、かなり大雑把にですが使い心地を試してみようと思います。ただ、最後までお読みになればお分かりになるかと思いますが色々つまずいたところもありますので、解決策をご存知の方はどしどしご指摘くださると有難いですm(_ _)m
tfp.stsモジュールでやってみる
簡単なデータを用いたチュートリアルは記事とNotebookの両方が公開されています。
なのですが、これだとほぼ単変量に近いデータなのでこの記事では他変量のデータでやってみようと思います。
サンプルデータとして用いるのは、以前RStanや{bsts}の実験材料としても用いたこのデータです。7日周期のseasonalityとトレンド（二階差分）、そして3つの「広告投下ボリューム」変数による回帰成分の3つの要素からなる非線形時系列データです。
それではやってみましょう。実際に使ったPythonコードの全体をGitHubにJuPyter Notebookとして上げてありますので、面倒という方はそちらをご覧ください。
%matplotlib inline import matplotlib as mpl from matplotlib import pylab as plt import matplotlib.dates as mdates import seaborn as sns import collections import numpy as np import pandas as pd import tensorflow as tf import tensorflow_probability as tfp from tensorflow_probability import distributions as tfd from tensorflow_probability import sts from sklearn.preprocessing import StandardScaler scaler = StandardScaler() # Data import d = pd.read_csv("https://raw.githubusercontent.com/ozt-ca/tjo.hatenablog.samples/master/r_samples/public_lib/jp/sample_marketing.csv") forecast_steps = 16 # Normalize exogenous variables d.iloc[:, 1:] = scaler.fit_transform(d.iloc[:, 1:]) d_train = d[:-forecast_steps] d_test = d[-forecast_steps:] # Plot the original data fig = plt.figure(figsize=(12, 6)) ax = fig.add_subplot(1, 1, 1) ax.plot(d_train.cv, lw=2, label="training data"); # Dedicated plotting function def plot_forecast(d, validation_steps, forecast_mean, forecast_scale, forecast_samples, title, x_locator=None, x_formatter=None): """Plot a forecast distribution against the 'true' time series.""" colors = sns.color_palette() c1, c2 = colors[0], colors[1] fig = plt.figure(figsize=(12, 6)) ax = fig.add_subplot(1, 1, 1) num_steps = len(d) num_steps_forecast = validation_steps num_steps_train = num_steps - num_steps_forecast ax.plot(d.iloc[:, 0], lw=2, color=c1, label='ground truth') forecast_steps = np.arange( num_steps_train, num_steps_train+num_steps_forecast) ax.plot(forecast_steps, forecast_samples.T, lw=1, color=c2, alpha=0.1) ax.plot(forecast_steps, forecast_mean, lw=2, ls='--', color=c2, label='forecast') ax.fill_between(forecast_steps, forecast_mean-2*forecast_scale, forecast_mean+2*forecast_scale, color=c2, alpha=0.2) ymin, ymax = min(np.min(forecast_samples), np.min(d.iloc[:, 0])), max(np.max(forecast_samples), np.max(d.iloc[:, 0])) yrange = ymax-ymin ax.set_ylim([ymin - yrange*0.1, ymax + yrange*0.1]) ax.set_title("{}".format(title)) ax.legend() if x_locator is not None: ax.xaxis.set_major_locator(x_locator) ax.xaxis.set_major_formatter(x_formatter) fig.autofmt_xdate() return fig, ax # Define model function like {bsts} def build_model(observed_time_series): trend = sts.LocalLinearTrend(observed_time_series=observed_time_series.iloc[:, 0]) seasonal = tfp.sts.Seasonal( num_seasons=7, num_steps_per_season=1, observed_time_series=observed_time_series.iloc[:, 0]) regression = sts.LinearRegression(design_matrix=observed_time_series.iloc[:, 1:]) model = sts.Sum([trend, seasonal, regression], observed_time_series=observed_time_series.iloc[:, 0]) return model # Estimate a model d_model = build_model(d_train) samples, kernel_results = tfp.sts.fit_with_hmc(d_model, d_train) # Extract MCMC samples with tf.Session() as sess: sess.run(tf.global_variables_initializer()) samples_, kernel_results_ = sess.run((samples, kernel_results)) # Forecast d_forecast_dist = tfp.sts.forecast(d_model, d, parameter_samples=samples_, num_steps_forecast=forecast_steps) d_forecast_mean = d_forecast_dist.mean()[..., 0] # shape: [50] d_forecast_scale = d_forecast_dist.stddev()[..., 0] # shape: [50] d_forecast_samples = d_forecast_dist.sample(10)[..., 0] # shape: [10, 50] # Extract forecast results with tf.Session() as sess: forecast_mean, forecast_scale, forecast_samples = sess.run(( d_forecast_mean, d_forecast_scale, d_forecast_samples)) # Extract regression coefficients for (param, param_draws) in zip(d_model.parameters, samples_): if param.name == "LinearRegression/_weights": weights = np.mean(param_draws, axis=0) print param.name, np.mean(param_draws, axis=0) # Compute regression component from the future values forecast_regression = np.dot(d_test.iloc[:, 1:], weights) # Plot the final result fig, ax = plot_forecast(d, forecast_steps, forecast_mean + forecast_regression, forecast_scale, forecast_samples, "Validation plot")
明らかにうまくいってませんorz そもそもこんなにマイナス方向に振れるはずがないんですが。。。一応MCMCサンプルのヒストグラムもチェックしたりしてみたものの、特に極端におかしい結果には見えない感じでした。
困ったこと：未来の説明変数を使った予測ができない
ところで、上のコードのこの箇所を見て「何だこれ？」と思った人もいるかもしれません。
# Compute regression component from the future values forecast_regression = np.dot(d_test.iloc[:, 1:], weights)
回帰成分の未来予測のところを別に計算しています。
というのも、tfp.sts.forecastの説明のどこを見ても未来の説明変数を与える方法が書いてないのです。それで僕は一旦回帰成分だけ別に計算する方法を取ったのですが、GitHubを見に行ったら同じポイントについて既にユーザーから指摘済みでissueも立っていて、開発チームの説明としては「回帰成分のdesign matrixを与えるところにtf.concat([...], axis=0)で未来の説明変数を与えろ」ということが書かれています。
その通りにやってみた結果がこちらなんですが、やっぱりうまくいっていませんorz
Rの{bsts}パッケージで同じことをやってみる
言うまでもなく、全く同じことがRの{bsts}パッケージでも出来ます。ということでやってみました。
d <- read.csv('https://raw.githubusercontent.com/ozt-ca/tjo.hatenablog.samples/master/r_samples/public_lib/jp/sample_marketing.csv') d_train <- d[1:84,] d_test <- d[85:100,] library(bsts) ss <- AddLocalLinearTrend(list(),d_train$cv) ss <- AddSeasonal(ss, d_train$cv, nseasons = 7) model <- bsts(cv~., state.specification = ss, niter=1000, data=d_train) # 略 pred <- predict(model, burn=100, newdata=d_test[,-1]) plot(pred, ylim = c(0, 1500), xlab = '', ylab = '') par(new = T) plot(d_test$cv, ylim = c(0,1500), xlim = c(-84,16), axes = F, type = 'l', col = 'red', lwd = 3, lty = 1, xlab = '', ylab = '') legend('topleft', legend = c('Predicted', 'Actual'), col = c('blue', 'red'), lty = 1, lwd = 3, cex = 1.5)
