状態空間モデルのれんしゅう（レストランの来客予測）

はじめに

時系列解析を行う上で、状態空間モデルが現状使用できる手法の中で優れているのではと自分の中で話題だったので、 使ってみることにしました。

状態空間モデルを使うことの目的

時系列解析を行う上でもモデルの種類は無数にあるのですが、 どういう場合に状態空間モデルを使うのでしょうか？

この部分については以下のサイトに詳しく書かれていました。
https://logics-of-blue.com/%E3%81%AA%E3%81%9C%E7%8A%B6%E6%85%8B%E7%A9%BA%E9%96%93%E3%83%A2%E3%83%87%E3%83%AB%E3%82%92%E4%BD%BF%E3%81%86%E3%81%AE%E3%81%8B/

状態空間モデルはデータが得られるプロセスとは何か？を重視したモデルであり、 回帰分析よりも各要因の影響が納得性の高い形で見ることができ、またその構造が 機械学習手法よりも理解しやすいものとなるようです。

説明性ではなく予測精度だけを追求したい場合は、統計モデルを使うことの優位性はないと思っています。（あればkaggleで使われているはずなので）
kaggleでは時系列系のコンペでも、過去データから特徴量を大量に作成し、勾配ブースティング系の学習器につっこむという方法が主流になっています。

状態空間モデルでは各説明変数がどう作用しているかを、理解できる形で出せることがメリットなので、 例えばレストランの売り上げについて考える場合、 曜日による影響や、複数ある広告媒体のそれぞれの影響等が把握できるようになり、 施策の打ち手を考えたり、商品開発の方向性を決めるときなんかに使える気がします。

使用データについて

今回はkaggleのRecruit Restaurant Visitor Forecastingのデータを使ってレストランの来客数の要因分解を行います。

このデータセットでは、 829店舗のレストランの2016/1/1-2017/04/22の期間の１日毎の売り上げが記載されています。 店舗には居酒屋やイタリアン等のジャンルの情報が入っており、 今回はジャンルが居酒屋である197店舗のみを対象とします。
（ジャンル毎にパラメータを設定したら収束しなくなったため）

事前調査

モデルを作る前にデータを分布を確認します。 下記は、日時に対する全レストランへの来客数の合計値の推移になります。2016/7で大きな値の変化がありますが、 これは登録されているレストランの数が増えていることが要因のようです。

次に、目的変数となる店舗毎の来客数の平均値のヒストグラムを確認します。

正規分布にするにはデータは少し左寄りですね。0以上ですし、ポアソン分布か対数正規分布で考えた方が良さそうです。

対数をとったグラフも出してみます。

少し微妙ですが、それほど正規分布とは外れていないと思いますので、 今回は目的変数の分布は対数正規分布としてモデルを組むことにします。

モデルの作成

stanのモデルを作成する上では下記のサイトを参考にさせていただきました。

• https://fisproject.jp/2017/10/state-space-model-with-stan/
• http://statmodeling.hatenablog.com/entry/state-space-model-many-terms
• https://logics-of-blue.com/multivariate-time-series-with-stan/

モデルの特徴量としては、トレンドと、週特徴、祝日特徴を考え、これらの特徴は店舗毎に共通であると仮定します。 トレンドの設定方法はいくつか方法がありますが、今回は前回の値のみを参照する形とします。 週特徴については参考のサイトのものをそのまま使わせてもらい、１週間の合計値が0になるようなパラメータの置き方をしました。 祝日特徴は、平日が祝日になるケースのときのみフラグを立てています。

このデータでは定休日の部分は欠損となっており、これをどう扱うかは悩みどころでした。 学習に使うのが1店舗だけであれば除外するという方法も考えられましたが、複数店舗で行う良いやり方が思いつかなかったため、 欠損の部分は１日前の値をそのまま使うことにしました。 （曜日効果の部分に影響が出る可能性がありますが、多店舗で実施するためその影響が小さくなってくれるのではという期待込みです。）

stanコードはこのようになりました。

model = """
    data {
        int n; // サンプルサイズ
        int week_day[n]; // 曜日
        int holiday_val[n]; //土日を除いた祝日フラグ
        int store_num; //店舗数
        matrix[n, store_num] y; //観測値
    }
    parameters {
        vector[6] s_raw;
        vector[store_num] r; //店舗毎のランダム効果
        real<lower=0> s_r;  //ランダム効果の標準偏差
        real muZero; // 左端
        vector[n] mu; // 確率的レベル
        real<lower=0> sigmaV; // 観測誤差の大きさ
        real<lower=0> sigmaW; // 過程誤差の大きさ
        real<lower=0> rh; //祝日効果
        
    }
    
    transformed parameters {
      vector[7] s;
      vector[n] week;
      
      s[1:6] = s_raw;
      s[7] = -sum(s_raw);
      for(i in 1:n){
         week[i] = s[week_day[i]+1];
         }
   
    }
    model {
        // 状態方程式
        mu[1] ~ normal(muZero, sqrt(sigmaW));
        for(i in 2:n) {
            mu[i] ~ normal(mu[i-1], sqrt(sigmaW));
        }
        
        // ランダム効果
        r ~ normal(0, s_r);
        
        // 観測方程式
        for(i in 1:n) {
            for(j in 1:store_num){
                y[i,j] ~ lognormal(mu[i] + week[i] + rh * holiday_val[i] + r[j], sqrt(sigmaV));
            }
        }
        
    }
"""
data={'n': visitors.shape[0], 'y': visitors.iloc[:,1:],'week_day':weeks,'holiday_val':holidays,
      'store_num':180}
fit = pystan.stan(model_code=model, data=data, iter=10000, chains=3)

結果

無事収束しました。

曜日毎の特徴を確認します。

# サンプル列を抽出
la  = fit.extract(permuted=True)
#曜日特性の確認
plt.figure(figsize=(12,8))
for i in range(la["s"].shape[1]):
    pd.Series(la["s"][:,i]).hist(label=i, alpha=0.5)
plt.legend()
plt.plot()

0-6が月-金に相当します。 金>=土>>>日>水>木>火>月という感じです。

次に、祝日効果を確認します。

#祝日効果の確認
plt.figure(figsize=(12,8))
pd.Series(la["rh"]).hist(bins=30)
plt.legend()
plt.plot()

0.04くらいが最頻のようです。 対数化しているので、実数に直すと1人相当です。 あまり影響が出ていないですね。。。

今回推定したパラメータの予測結果がどれだけ過去データにどれだけあてはまっているかをいくつか見てみます。

店舗サンプル１

pred=[]
for i in range(visitors.shape[0]):
    pred.append(np.exp(la["mu"][:,i].mean(axis=0)+la["week"][:,i].mean(axis=0)+la["rh"].mean()*holidays.values[i]+la["r"][:,0].mean()))
plt.figure(figsize=(12,8))
plt.plot(pred[:100])
plt.plot(visitors["air_0241aa3964b7f861"][:100])

ベースはあっていますが波形をほとんど追えていません。 今回全店舗で週特徴の効果が同様であるという仮定をおいているため、 この店舗のように週効果が他の店舗と異なるものには対応できなかったようです。

店舗サンプル２

pred=[]
for i in range(visitors.shape[0]):
    pred.append(np.exp(la["mu"][:,i].mean(axis=0)+la["week"][:,i].mean(axis=0)+la["rh"].mean()*holidays.values[i]+la["r"][:,1].mean()))
plt.figure(figsize=(12,8))
plt.plot(pred[:100])
plt.plot(visitors["air_03963426c9312048"][:100])

上下幅は小さいですが、先ほどの店舗と比較すると上下の変動を追従できている気がします。

まとめ

今回はあまり当てはまりの良いモデルとはなりませんでした。 全店舗で週効果が同じだと仮定したことがよくなかったかもしれませんし、 欠損値の補完方法についても検討の余地があると思います。 また、祝日効果がかなり低めに出てしまっていましたが、 日曜より土曜の方が来客が増えることを考えると祝日ではなく祝前日の設定をしておけばもう少しあてはまりがよくなったかなと思います。

メモ

• 各説明変数をどうモデルに入れるか、という部分は試行錯誤が必要。
• 試行錯誤が必要だが、stanは時間がかかる。記事のモデルだと２時間くらいかかっているが、lightgbmだったら同じくらいのデータ量でも数分もかからないと思うので、データ量が一定以上あるときは厳しさを感じる。
• 自分の考えた最強の状態空間モデルは大抵収束しない。

使用したコードは以下
https://github.com/rmizuta3/ssmodel_recruit