読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

[dlm]状態空間モデルでトレンドと広告の効果を分離して推定する

はじめに

Stanで統計モデリングを学ぶ(7): 時系列の「トレンド」を目視ではなくきちんと統計的に推定する - 東京で働くデータサイエンティストのブログ をみてください。

上記事では Stan で状態空間モデルを推定しているので、ここでは R の dlm パッケージを使ってやってみます。

モデル化

コンバージョン(CV)がベースラインとなるトレンドと三種類の広告の効果の和で表わされるというモデルです。

 \rm{CV}_{t}=\mu_t + \beta_1 {\rm ad}_{1,t}+ \beta_2 {\rm ad}_{2,t}+ \beta_3 {\rm ad}_{3,t}+d + v _t

ここで d は切片、v _t は誤差項で独立に正規分布に従います。

トレンド  \mu_t は二回差分で表わされ、

 \mu_t =2\mu _{t-1} - \mu_{t-2} + \eta_{1,t}

です。\eta _{1,t} は誤差項で独立に正規分布に従います。

この式の意味は、「今回の変化量と前回の変化量はだいたいおなじくらいだろう」という仮定です。

\mu _ t - \mu_{t-1} \simeq \mu _ {t-1} - \mu_{t-2}

\mu _ {t-1} を移行してやると、

\mu _ t\simeq 2\mu _ {t-1} - \mu_{t-2}

という形がでてきます。

行列表記

dlm が想定しているモデルは

y_t =F_t \theta_t +v_t, v_t \sim N(0,Vt)
\theta _t=G_t \theta_{t−1}+w_t, w_t\sim N(0,Wt)

という形なので上で仮定したモデルを行列形式に書きなおします。

\beta_1\beta_2\beta_3、それに d_t は時間に依存しないので、

\beta_{1,t} = \beta_{1,t-1}
\beta_{2,t} = \beta_{2,t-1}
\beta_{3,t} = \beta_{3,t-1}
d_{t} = d_{t-1}

と改めて表します。

 \rm{CV}_{t}=\begin{pmatrix}1&0 &\rm{ad}_{1,t}&\rm{ad}_{2,t}&\rm{ad}_{3,t}&1\end{pmatrix}
\begin{pmatrix}\mu_t \\ \mu_{t-1} \\ \beta_{1,t} \\  \beta_{2,t} \\ \beta_{3,t} \\ d_t \end{pmatrix}+v_t

 \begin{pmatrix}\mu_t \\ \mu_{t-1} \\ \beta_{1,t} \\  \beta_{2,t} \\ \beta_{3,t} \\ d_t \end{pmatrix}=\begin{pmatrix}2 & -1&0&0&0&0\\
                  1&0&0&0&0&0\\
                  0&0&1&0&0&0\\
                  0&0&0&1&0&0\\
                  0&0&0&0&1&0\\
                  0&0&0&0&0 &1\end{pmatrix}
\begin{pmatrix}\mu_{t-1} \\ \mu_{t-2} \\ \beta_{1,t-1} \\  \beta_{2,t-1} \\ \beta_{3,t-1} \\ d_{t-1} \end{pmatrix}+w_t

\beta _t d_t には誤差項がつきませんでしたが、これは「平均 0 分散 0 の正規分布に従う誤差がついている」と捉えます。

 W_t = \begin{pmatrix}
\eta_1&0&0&0&0&0\\
0&0&0&0&0&0\\
0&0&0&0&0&0\\
0&0&0&0&0&0\\
0&0&0&0&0&0
\end{pmatrix}

です。

dlm

上記のモデルを愚直に打ち込みます。

build1 <- function(theta){
  dlm(FF=matrix(c(1,0,1,1,1,1),1,6),
      V=matrix(exp(theta[1])),
      GG=matrix(c(2,-1,0,0,0,0,
                  1,0,0,0,0,0,
                  0,0,1,0,0,0,
                  0,0,0,1,0,0,
                  0,0,0,0,1,0,
                  0,0,0,0,0,1),6,6,byrow = TRUE),
      W=matrix(c(exp(theta[2]),0,0,0,0,0,
                 0,0,0,0,0,0,
                 0,0,0,0,0,0,
                 0,0,0,0,0,0,
                 0,0,0,0,0,0),6,6,byrow = TRUE),
      JFF=matrix(c(0,0,1,2,3,0),1,6),
      X = as.matrix(data1[,2:4]),
      m0 = c(0,0,0,0,0,0),
      C0 = matrix(c(1e+7,0,0,0,0,0,
                    0,0,0,0,0,0,
                    0,0,1e+7,0,0,0,
                    0,0,0,1e+7,0,0,
                    0,0,0,0,1e+7,0,
                    0,0,0,0,0,1e+7),6,6,byrow = TRUE))
}

FF が F_t、GG が G_t に対応します。
m0 は \theta_t の事前分布の平均、C0 が分散です。
JFF は FF のどの成分に説明変数が入るのかを指定しています。

data1 <-read.table("https://raw.githubusercontent.com/ozt-ca/tjo.hatenablog.samples/master/r_samples/public_lib/jp/stan_trend_issue.txt",header = TRUE)
fit1 <- dlmMLE(data1$cv, parm=c(log(var(data1$cv)),1), build1, method = "Nelder-Mead")

model1 <- build1(fit1$par)
Filt1 <- dlmFilter(data1$cv,model1)
Smooth1 <-dlmSmooth(Filt1)

plot(data1$cv,type="b")
lines(dropFirst(Smooth1$s[,1]) + apply(Smooth1$s[-1,3:5]*data1[,2:4],1,sum) + dropFirst(Smooth1$s[,6]),col="red",lwd=2)
lines(dropFirst(Smooth1$s[,1]+Smooth1$s[,6]),col="blue",lwd=2)

パラメータを推定して、スムージングしてやると下図のようにトレンドと広告の効果を分離することができました。

青い線がトレンド、赤い線がトレンドと広告の効果を合わせた CV の推定値です。

f:id:abrahamcow:20160507070621p:plain

trend <- read.table("https://raw.githubusercontent.com/ozt-ca/tjo.hatenablog.samples/master/r_samples/public_lib/jp/trend_actual.txt")
plot(cumsum(trend$V1),type="b")
lines(dropFirst(Smooth1$s[,1]),col="blue",lwd=2)

トレンドの正解データと推定されたトレンドをあわせてプロットしてやります。

f:id:abrahamcow:20160507071034p:plain

まずまず正解に近い値が求まっていると言えそうです。

参考文献

本エントリはかなり説明をはしょっているので、以下を参照していただけると幸いです。

状態空間時系列分析入門

状態空間時系列分析入門

R で 状態空間モデル: {dlm} で単変量モデルの状態空間表現を利用する - StatsFragments

dlmによる時変係数モデル | 状態空間モデル | Logics of Blue