カルマンフィルタのアドテクへの応用（実践編）
実験はR言語で行いました。ソース及びデータの一式はChiral's Gistに置いてあります。
具体的な問題設定
ある市場において、商品カテゴリ１，２があり、A社とB社が競合してるとします。A社は１，２両方のカテゴリの商品を扱っているいっぽうで、B社はカテゴリ２の商品のみを扱っています。
ユーザの興味ベクトルは、商品カテゴリ１，２それぞれの選好度、A社とB社それぞれのブランド選好度という４つの要素と、それらの4要素の変化の速度成分を加えて計8次元のベクトルとします。ここでいうユーザは何らかのひとかたまりのオーディエンスという想定です。
そしてA,Bの2社が、ユーザに向けてそれぞれ広告を打ちました。時間の単位は1週間とし、両社の広告の応酬が繰り広げられた3か月(13週)間を分析対象とします。
広告の種類は3種類、すなわち広告１がA社の商品カテゴリ１、広告２がA社の商品カテゴリ２、広告３がB社の商品カテゴリ２、です。
広告配信量は一週間ごとに１０００インプレッションとします。1週間ごとに広告クリック数を集計したものを観測データとします。
さてここで、広告配信について次のようなシナリオを想定してます。
- 最初の5週間は、A社が商品カテゴリ１（広告１）を広告配信した。
- 続く4週間で、B社が商品カテゴリ２（広告３）を広告配信した。
- 最後の4週間は、またA社が商品カテゴリ２（広告２）を広告配信した。
クリック数のデータも合わせて手入力でそれっぽく作ったデータが以下の図になります。左がインプレッション、右がクリック数、横軸はいずれも週単位の時間軸です。
黒線がA社、赤線がB社の広告及び反応です。A社は計9週間で９０００インプレッション、B社は4週間で４０００インプレッションと、使った広告予算には2倍以上の開きがあるという問題設定であることに注意しておきます。
Rプログラム
具体的な行列のパラメータ設定や計算手順を述べると長くなりそうなのでRのソースを貼り付けておきます。結果のグラフと考察は次節にて。
# 単位時間間隔Δt
dt <- 1
# 興味ベクトルの仕様
# カテゴリ選好度2つとブランド選好度2つの4次元と
# それらの速度成分を加えた合計8次元
# c(cat1,cat2,brandA,brandB,velo1,velo2,veloA,veloB)
interest_elems <- list("cat1","cat2","brandA","brandB")
dimInterestHalf <- length(interest_elems)
dimInterest <- dimInterestHalf * 2
# バナーは3種類とする
# 各バナーに興味ベクトルを付与
dimAd <- 3
# スケールを維持するため合計が1になるように分配
ad1 <- c(1/2,0,1/2,0) # ad1 = cat1 + brandA
ad2 <- c(0,1/2,1/2,0) # ad2 = cat2 + brandA
ad3 <- c(0,1/2,0,1/2) # ad3 = cat2 + brandB
# 行列H： 興味ベクトルxから反応ベクトルzへの変換行列
velo0 <- rep(0,dimInterestHalf)
H <- t(matrix(c(ad1,velo0,
ad2,velo0,
ad3,velo0),
ncol=dimAd, nrow=dimInterest))
# 行列R： 観測ベクトルzの観測ノイズの共分散行列（対称行列）
R <- diag(0.5 * c(1,1,1))
# 行列F1： 興味ベクトルxの時間遷移行列
F1_00 <- diag(dimInterestHalf)
F1_10 <- matrix(0,dimInterestHalf,dimInterestHalf)
F1_01 <- dt * diag(dimInterestHalf)
F1_11 <- diag(dimInterestHalf)
F1 <- cbind(rbind(F1_00,F1_10),rbind(F1_01,F1_11))
# 行列F2： 広告ベクトルuから興味ベクトルxへの変換行列
alpha <- 0.01
F2 <- alpha * matrix(c(ad1,ad2,ad3),ncol=dimAd, nrow=dimInterestHalf)
# 行列Q： 興味ベクトルの単位時間あたりのノイズ源 w_t
Q <- diag(0.025 * c(1,1,1,1))
# 行列F3： ノイズ源w_t を興味ベクトルへの寄与に変換する行列
F3_0 <- (1/2*dt*dt) * diag(dimInterestHalf)
F3_1 <- dt * diag(dimInterestHalf)
F3 <- rbind(F3_0,F3_1)
F2 <- F3 %*% F2 # 広告インプレッションを力の作用とみなすため
# ノイズの寄与と同じ変換を掛ける
# 初期条件 x0,P0
x0 <- rep(0,dimInterest)
P0 <- diag(dimInterest)
x <- rbind(x0)
P <- list(P0)
# データの読み込み
z <- read.csv("response.csv")
u <- read.csv("impression.csv")
# 予測ステップ
# 注）R言語の添え字の都合により、、
# x,P と z,u は時刻インデックスが一つずれてる
kalman_predict <- function(t) {
x1 <- F1 %*% x[t,] + F2 %*% t(u[t,])
x <<- rbind(x,t(x1))
P[[t+1]] <<- F1 %*% P[[t]] %*% t(F1) + F3 %*% Q %*% t(F3)
}
# 更新ステップ
kalman_update <- function(t) {
e <- z[t,] - x[t+1,] %*% t(H)
S <- R + H %*% P[[t+1]] %*% t(H)
K <- P[[t+1]] %*% t(H) %*% solve(S)
I <- diag(dimInterest)
x[t+1,] <<- x[t+1,] + K %*% t(e)
P[[t+1]] <<- (I - K %*% H) %*% P[[t+1]]
}
# シミュレーション
upper <- c()
lower <- c()
for (t in 1:nrow(z)) {
print(paste("time =",t))
kalman_predict(t)
print(paste(" predict_x =",x[t,]))
kalman_update(t)
print(paste(" update_x =",x[t,]))
sd <- sqrt(diag(P[[t]]))
upper <- rbind(upper,x[t,]+sd)
lower <- rbind(lower,x[t,]-sd)
}
# 結果のプロット
par(mfrow=c(3,2))
xlim<-c(1,nrow(z))
ylim<-c(min(u),max(u))
plot(u[,1],type='l',col='1',ylab='u',xlim=xlim,ylim=ylim)
par(new=T)
plot(u[,2],type='l',col='1',ylab='',xlim=xlim,ylim=ylim)
par(new=T)
plot(u[,3],type='l',col='2',ylab='',xlim=xlim,ylim=ylim)
ylim<-c(min(z),max(z))
plot(z[,1],type='l',col='1',ylab='z',xlim=xlim,ylim=ylim)
par(new=T)
plot(z[,2],type='l',col='1',ylab='',xlim=xlim,ylim=ylim)
par(new=T)
plot(z[,3],type='l',col='2',ylab='',xlim=xlim,ylim=ylim)
for (i in 1:dimInterestHalf) {
ylim<-c(min(lower[,i]),max(upper[,i]))
plot(x[,i],type='l',ylab=interest_elems[[i]],xlim=xlim,ylim=ylim)
par(new=T)
plot(upper[,i],type='l',ylab='',col='4',xlim=xlim,ylim=ylim)
par(new=T)
plot(lower[,i],type='l',ylab='',col='2',xlim=xlim,ylim=ylim)
}
結果
上のプログラムに13週間ぶんのデータを食わせると、以下のように興味ベクトルの要素ごとに時間推移をプロットしてくれます。
黒い線が各時刻での更新後の値、青と赤の線は誤差行列
の対角成分をそれぞれ平方根を取って標準偏差的なスケールにして
の上下の範囲を計算したものです。縦軸は、クリックレートの1000倍のスケールになるようにパラメータ設定しています。例えば、２５だったらCTR=2.5%といった意味合いになります。
理論編でも述べましたが、モデルの中でガウスノイズの数式は出てくるのですが、プログラムで計算することは乱数なしの決定的な計算になるのがカルマンフィルタ。なので乱数頼みのモンテカルロ法や初期値の影響を受けるEM学習等とちがって、推定結果が解析的にビシッと一つに求まるのでやってみて実にすがすがしいです。パラメータを色々与えないといけないところが難点ですがそれもこの結果のグラフで報われる感じもします。
さて、グラフの解釈をしましょう。グラフの上２つのは商品カテゴリ１，２の選好度の推移、下２つはブランドA,Bの選好度の推移です。左↑のカテゴリ１はまぁ普通ですが右のカテゴリ２が途中6週目あたりで不思議にもマイナスになっています。これはつまり、最初の5週間は両方のカテゴリの商品を持つA社がカテゴリ１のみ広告を打っておりそのおかげでA社のブランド選好度が向上したのでそのままであればA社が扱っているカテゴリ２の反応も向上するはずが実際の観測データではゼロだったのでつり合いを取るためにカテゴリ２自体の選好度をマイナスに持っていった、ということだと思われます。
そして面白いことにB社のほうは最初の5週間でまったく広告を打っていないにも関わらずその間のブランド選好度は向上しています。カテゴリ２の選好度が下がったにも関わらずB社のカテゴリ２への反応データがゼロだったのでつり合いを取るためにB社のブランドをプラスに持っていったということだと思われます。
このあたり、ブランドと商品カテゴリの各成分への分解をカルマンフィルターが自動的にやってくれていて非常に面白いです。
続く6週目で今度はB社が広告配信を開始しますが、6週目から8週目までのB社のブランドはかなり効率よく向上していってます。最初の5週間でA社が広告を打ったときは０→２０くらいだったのに、続く4週間のB社の広告は１０→５０も向上しています。これは観測データを恣意的に設定したのではないことを確認するため、上のほうで掲載した反応データの図を再掲します。
むしろB社の広告（赤線）への反応は総じて小規模になってます。
ではなぜB社のブランド力は向上したのか。それは、まず6週目の時点で商品カテゴリ２の選好度はマイナスであり、A社の最初の広告配信は終わったのでA社からのカテゴリ２への寄与はこれ以降減少の一途です。するとB社のカテゴリ２の広告としては、カテゴリ2自体の選好度でもって反応を期待することはできなくなり、B社のブランドとしての訴求力しか広告への反応を支えるものがありません。そして実際に広告の反応はあって配信中は反応が向上しました。そのことを持ってカルマンフィルターはこれは殆どがB社のブランド力の向上の賜物であると推定したというわけです。
これも非常に面白い推定結果です。カテゴリ自体の認知度が低く、市場を盛り上げてくれる競合もいないときは、たとえ小さくても市場からの反応があればブランド資産は蓄積されていってるのだということを定量的にモデル化できたことを示唆しているようにも思います。
さて、つづいて最後の4週間ですがB社の広告配信が終わって今度はA社がカテゴリ２の広告を打ちました。カテゴリ２の認知・選好度が向上してA社のブランド選好度が下げ止まって上昇に転じているという、ここは大体予想通りの結果です。
最終的に、A社は９０００インプレッション、B社は４０００インプレッションという広告予算投下において、観測データを見るとA社はコスト掛けただけあって良い反応が得られていますが、ブランド選好度を見ると乱高下して場当たり的な感もあります。刈取りをしている感じでしょうか。B社のほうはカテゴリ２の選好度が低いときに広告を打ったので反応はイマイチでしたがブランド資産の蓄積がグラフに表れました。もしこのカルマンフィルターのモデルが正しければブランド資産を定量指標としてその後のマーケティングプランニングの指針になることでしょう。
感想＆所感
wikipediaをサーフィンしながら着想を思いついたのが昨日の午後で、（途中遊びにいって帰ってから徹夜してやりましたが）理論の適用からRプログラミングと実験、ブログ執筆まで含めて18時間程度でやったのでなかなかやりごたえありました。結果も少し示唆を感じるものが出てきたので面白かったです。
ただ、カルマンフィルタは状態遷移と観測のモデル（行列）のパラメータ数が多く、またモデルの正確さに依拠して解析的にバチンと答えを出す手法なので、ちょっとカタい感じがしました。モデルが正確ならよいですが、パラメータ設定の恣意性は免れません。
それと、今回は対象オーディエンスを集団としましたが、カルマンフィルタで一人一人の興味ベクトルを計算しようとすると頻度の少ない確率事象を線形モデルで表現しずらいという問題があることがわかったので、非線形な粒子フィルタを引き続き検討しようと思いました。
カルマンフィルタの使い方は結構分かったので、他の課題に適用できそうなときがあったらまたやってみたいと思います。
