廿TT

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

ゼロの多いコンバージョンの状態空間モデル

今日の川柳

対象

アマゾンアフィリエイトをやってる。

コンバージョンはあまりない。

いくつか図を貼ります。

去年一年間の日ごとの注文数合計です。

f:id:abrahamcow:20180613033935p:plain

なにもない日がけっこう多い。でも後半ちょっと黒いところが増えているような気もする。

と思ったけど、以下を見るとやっぱりコンバージョン数は年間通じて一定かな。

f:id:abrahamcow:20180613035606p:plain

コンバージョンのレートはどうだろう。

f:id:abrahamcow:20180613035629p:plain

なにもない日が多くてスパイク状に突き出して見える。

時間の情報を潰してヒストグラムにする。

f:id:abrahamcow:20180613035722p:plain

成功確率一定の二項分布にしては0が多い。また幅が広い。

そこで以降でコンバージョンレートのトレンドを推定して見ることにします。

ちないに、図を書いたRのコードは以下です。

library(tidyverse)
library(tcltk)

dat <-read_csv("~/Downloads/report2017.csv",skip = 1)

theme_set(theme_classic(base_size = 14, base_family = "Osaka"))

ggplot(dat,aes(x=日付,ymin=0,ymax=注文数合計))+
  geom_linerange()

ggplot(dat,aes(x=日付,y=cumsum(注文数合計)))+
  geom_step()

ggplot(dat,aes(x=日付,y=注文数合計/クリック数))+
  geom_line()

ggplot(dat,aes(x=注文数合計/クリック数))+
  geom_histogram(bins=20,fill="white",colour="black")

データは
report2017.csv · GitHub
に置いておきます。

モデル

以下のモデルを考えます。

Yang (2012) http://ir.uiowa.edu/cgi/viewcontent.cgi?article=3166&context=etd では zero-inflated ポアソン分布を使った状態空間モデルを論じています。

この論文の内容はRのパッケージとしても提供されています(CRAN - Package ZIM)。

これを参考に二項分布のモデルを考えます。 y_t がコンバージョン数、n_t がクリック数です。ほかは潜在変数です。ここで t=1,\ldots,T とします。

 y_t \mid z_t=1 \sim {\rm binomial}(n_t, {\rm logit}^{-1}(s_t))
 y_t \mid z_t=0 \sim \mbox{0 with probability 1}
 z_t \sim {\rm binomial}(1, \phi)
 s_t \mid s_{t-1} \sim {\rm normal}(s_{t-1}, \sigma^2)
 s_{0} \sim {\rm normal}(0, 1)

最初の時点は情報がないので無情報事前分布で「ピンどめ」してやる必要があります。

s はロジットスケールで動くので、標準正規分布でも十分幅が広いと考えました。

標準正規乱数を逆ロジット変換するとこんな感じになります。

f:id:abrahamcow:20180613040046p:plain

hist(plogis(rnorm(10000,0,1)),breaks = "scott")

0.1とか0.9とか、高い確率や低い確率もそこそこカバーされています。

パーティクルスムーザー

フィルタリング分布とは t 時点までの情報を使った t 時点の状態変数の分布  p(s_{t}|y_{1:t}) を指します。

スムージング(平滑化)分布とは最後の T 時点までの情報を使った t 時点の状態変数の分布  p(s_{t}|y_{1:T}) を指します。

フィルタリングにはパーティクルフィルタを使います。

パーティクルフィルタのアルゴリズムは、R によるすごくかんたんなパーティクルフィルタの実装例 - 廿TT に書いたことがある。

  1. 一期先の予測を乱数でばらまく
  2. 得られた乱数に尤度の重みをつけてまき直す

の繰り返しです。

スムージングには Godsill et al. (2004) http://ftp.stat.duke.edu/WorkingPapers/03-12.pdf の方法を使います。

マルコフ性の仮定より、

 p(s_t|s_{t+1:T}, y_{1:T}) = p(s_{t}|s_{t+1}, y_{1:T})

この確率密度はベイズの定理より以下に比例します。

 p(s_{t}|y_{1:t})p(s_{t+1}|s_{t})

一つ目の因子は尤度、フィルタリングの際に使う重みと一緒です。二つ目は状態変数の尤度みたいなものです。

これを重みにしてパーティクルをリサンプリングすることで平滑化分布が得られます。

モンテカルロEMアルゴリズム

パーティクルスムーザーで未観測の変数が求まったらその標本平均を使ってパラメータを更新していけばよい。

二項分布のパラメータの最尤推定量は標本平均、正規分布の分散パラメータの最尤推定量は標本分散です。

R のコード

PF <- function(y,n,sigma,phi,np=10000){
  s <- matrix(0,366,np)
  z <- matrix(0,365,np)
  wt <- wt2 <-matrix(0,365,np)
  s[1,] <- rnorm(np,0,1)
  ll <- 0
  for(i in 1:365){
    s[i+1,] <- rnorm(np,s[i,],sigma)
    z[i,] <-rbinom(np,1,phi)
    wt[i,] <-dbinom(y[i],n[i],z[i,]*plogis(s[i+1,]),log = TRUE)
    ind <-sample.int(np, replace = TRUE, prob = softmax(wt[i,]))
    s[i+1,] <- s[i+1,ind]
    z[i,] <- z[i,ind]
    ll <- ll + logmeanexp(wt[i,])
  }
  for(i in 365:1){
    wt2[i,] <- logsumexp(wt[i,],dnorm(s[i+1,],s[i,],sigma,log=TRUE)+
                           dbinom(z[i,],1,phi,log=TRUE))
    ind <-sample.int(np, replace = TRUE, prob = softmax(wt2[i,]))
    s[i,] <- s[i,ind]
    z[i,] <- z[i,ind]
  }
  return(list(s=s,z=z,ll=ll))
}

softmax、logsumexp、logmeanexpはオーバーフローしないようにちょっと工夫している。

softmax <- function(x){
  tmp <- max(x,na.rm = TRUE)
  out <-exp(x-tmp)/sum(exp(x-tmp),na.rm = TRUE)
  out[is.na(out)] <- 0
  return(out)
}
logsumexp <- function(x,y){
  x + log(1+exp(y-x))
}
logmeanexp <- function(x){
  tmp <- max(x)
  tmp + log(sum(exp(x-tmp))) - log(length(x))
}

EMアルゴリズムの部分は以下のようにしました。

R の sd は本当は最尤推定量ではないけど、まあたぶん影響ないと思うのでそのままです。

phi <- 0.1
sigma <-0.1
llhist <- numeric(10)
pb <- txtProgressBar(min = 1, max = 10, style = 3)
set.seed(1234)
for(k in 1:10){
  filt1 <-PF(y = dat$注文数合計, n = dat$クリック数, phi = phi, sigma = sigma)
  phi <- mean(filt1$z)
  sigma <-sd(diff(rowMeans(filt1$s)))
  llhist[k] <- filt1$ll
  setTxtProgressBar(pb, k)
}

10回しかまわしてないけど尤度は十分大きくなっていると判断しました。

f:id:abrahamcow:20180613043559p:plain

ggplot(NULL,aes(x=1:10,y=llhist))+
  geom_point()+
  geom_line()+
  labs(x = "iter", y = "log-lik")

結果

コンバージョンレートのトレンドは以下のようになりました。

f:id:abrahamcow:20180613043659p:plain

周期性はもっと一週間単位とか短いかと思っていたけど一ヶ月くらい? で動いている。

長期でみると上昇傾向が見られます。

phat <-plogis(rowMeans(filt1$s)[-1])
CI <-plogis(apply(filt1$s[-1,],1,quantile,prob=c(0.025,0.975)))

ggplot(dat,aes(x=日付,y=注文数合計/クリック数))+
  geom_line()+
  geom_line(aes(y=phat),colour="royalblue",size=1)+
  geom_ribbon(aes(ymin=CI[1,],ymax=CI[2,]),fill="royalblue", alpha=0.3)

φ と σ の推定値は以下です。

> print(phi)
[1] 0.2843312
> print(sigma)
[1] 0.0343847

トレンドは上昇傾向でもコンバージョンがあるのは3日か4日に一回くらいみたいです。

感想

パーティクルスムーザーで状態推定してEMアルゴリズムでパラメータ推定というテクニックはけっこう便利そう。

Rcpp とか使えば速くなるだろうし実装もそんなに難しくない(場合がある)。

アマゾンアフィリエイトのコーナー

最近読んでおもしろかった本です。

あと、この文章がわかりにくいと感じた方には以下の本がおすすめです。