六本木で働くデータサイエンティストのブログ

元祖「銀座で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木

TensorFlow Probabilityのtfp.stsモジュールを使って構造時系列モデリングを回してみる

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")

f:id:TJO:20190427123622p:plain

明らかにうまくいってません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)

f:id:TJO:20190427163349p:plain

こちらではうまくいきました。ということでtfp.stsについては要経過観察といったところでしょうか。。。