Restricted Boltzmann Machineを実装してみた
R言語で実装してみました。特にパッケージは使ってないのでコピペすれば動くと思います。
このプログラムは隠れノード数やCD-kのkや更新ステップのηや収束条件のパラメータ等を変えて試すことが出来ます。
また学習後のRBMを使ってreconstruct（後述）したり隠れノードのサンプリングや確率計算をしたり出来ます。
ソース
### Restricted Boltzmann Machine implementation by isobe
sigmoid <- function(x) 1/(1+exp(-x))
rbm <- function(obs,n_hidden,eta=0.05,
epsilon=0.05,maxiter=100,
CD_k=1,reconstruct_trial=10,
verbose=0) {
L <- nrow(obs)
N <- ncol(obs)
M <- n_hidden
# initial values assinment
# cf) Chapter 8 in
# http://www.cs.toronto.edu/~hinton/absps/guideTR.pdf
pn <- apply(obs,2,function(x) min(0.9,sum(x)/L))
bn <- log(pn/(1-pn))
bm <- rep(0,M)
W <- matrix(rnorm(N*M,0,0.01),N,M)
pv_h <- function(i,h) {
sigmoid(sum(W[i,]*h)+bn[i])
}
ph_v <- function(i,v) {
sigmoid(sum(W[,i]*v)+bm[i])
}
gs_step <- function(x,n,p_func) {
r<-c()
for (i in 1:n) {
r<-c(r,rbinom(1,1,p_func(i,x)))
}
return(r)
}
gs_v <- function(h) gs_step(h,N,pv_h)
gs_h <- function(v) gs_step(v,M,ph_v)
cd_k <- function(v) {
v1 <- v
for (i in 1:CD_k) {
h1 <- gs_h(v1)
v1 <- gs_v(h1)
}
# R has immutable value and lexical scope,
# so we can overwrite locally.
for (i in 1:N) for (j in 1:M) {
W[i,j] <- ph_v(j,v)*v[i]-ph_v(j,v1)*v1[i]
}
bn <- v-v1
for (j in 1:M) {
bm[j] <- ph_v(j,v)-ph_v(j,v1)
}
return(list(W=W,bn=bn,bm=bm))
}
theta_step <- function() {
W <- matrix(0,N,M)
bn <- rep(0,N)
bm <- rep(0,M)
for (i in 1:L) {
if (verbose>=3) cat(paste("theta for obs ",i,"\n"))
d <- cd_k(obs[i,])
W <- W+d$W
bn <- bn+d$bn
bm <- bm+d$bm
}
return(list(W=W,bn=bn,bm=bm))
}
reconstruct <- function(v) gs_v(gs_h(v))
recon_error <- function() {
r <- 0
for (t in 1:reconstruct_trial) for (i in 1:L) {
v <- obs[i,]
v1 <- reconstruct(v)
r <- r+sum(abs(v-v1))
}
return(r/(N*L*reconstruct_trial))
}
err <- 1
count <- 0
cat("init OK. \n")
while (err>epsilon && count<maxiter) {
if (verbose>=2) cat(paste("step =",count,"\n"))
d <- theta_step()
backup <- list(W=W,bn=bn,bm=bm,err=err)
W <- W + eta*d$W
bn <- bn + eta*d$bn
bm <- bm + eta*d$bm
count <- count+1
err <- recon_error()
if (backup$err<err) {
W <- backup$W
bn <- backup$bn
bm <- backup$bm
err <- backup$err
} else if (verbose) {
if (verbose>=1) print(paste("step",count,": err=",err))
}
}
hidden_prob <- function(v) {
apply(rbind(1:M),1,function(i) ph_v(i,v))
}
learn_info=paste("step",count,": err=",err)
obj <- list(W=W,bn=bn,bm=bm,
learn_info=learn_info,
hidden_prob=hidden_prob,
hidden_sample=gs_h,
reconstruct=reconstruct)
class(obj) <- 'rbm'
return(obj)
}
print.rbm <- function(rbm) {
cat("edge weights:\n")
print(rbm$W)
cat("\nbias for observable nodes:\n")
print(rbm$bn)
cat("\nbias for hidden nodes:\n")
print(rbm$bm)
cat(paste("\n",rbm$learn_info,"\n",sep=''))
}
rbm_hidden_prob <- function(obj,obs) obj$hidden_prob(obs)
rbm_hidden_sample <- function(obj,obs) obj$hidden_sample(obs)
rbm_reconstruct <- function(obj,obs) obj$reconstruct(obs)
### test program
test <- function() {
obs <- rbind(c(1,0,1),
c(1,1,0),
c(1,0,1),
c(0,1,1))
net <- rbm(obs,2,verbose=1,maxiter=3000)
print(net)
x <- c(1,1,0)
trial <- 5
cat("original")
print(x)
for (t in 1:trial) {
cat("reconstructed")
print(rbm_reconstruct(net,x))
}
}
test()
テスト
上のソースには観測ノード数が3つ、隠れノード数が3つでXORすると1になるパリティチェック的なものを構成できるかのテストが書いてあります。実行すると以下のような感じに。
init OK.
[1] "step 1 : err= 0.408333333333333"
[1] "step 2 : err= 0.35"
[1] "step 46 : err= 0.35"
[1] "step 53 : err= 0.35"
[1] "step 104 : err= 0.341666666666667"
[1] "step 105 : err= 0.325"
[1] "step 158 : err= 0.325"
[1] "step 321 : err= 0.325"
[1] "step 363 : err= 0.283333333333333"
[1] "step 444 : err= 0.283333333333333"
[1] "step 1181 : err= 0.266666666666667"
edge weights:
[,1] [,2]
[1,] 0.07643391 0.08229904
[2,] 0.02417799 0.02915906
[3,] -0.05989687 -0.07663989
bias for observable nodes:
[1] 1.2486123 0.0500000 0.9486123
bias for hidden nodes:
[1] -9.421639e-04 -5.015886e-05
step 3000 : err= 0.266666666666667
original[1] 1 1 0
reconstructed[1] 0 0 1
reconstructed[1] 1 1 1
reconstructed[1] 1 0 1
reconstructed[1] 1 0 1
reconstructed[1] 1 1 13ビットのパターンでパリティが0になる4つを与えて、学習させた結果で110というビットを食わせて観測ノード→隠れノード→観測ノードとしたときの確率でサンプリング（reconstructionというそうです）を何度かやってみるというものです。
RBMはそもそも1層で何かを学習させるものではないと思うのでまぁあんまりうまくいってないですね。他の関数として8bitの2進数表現で4以上の数値を与えるようにするとパリティ関数よりは少し精度があがりました。
というわけで次はこれを使ってDeep Belief Networkにチャレンジしてみます。ワクワク！
ロジットとプロビットの使い分け
出力変数（被説明変数）がYes/Noみたいな２値で表されるようなモデルを学習させたい場合についてググるとロジスティック回帰とかプロビット回帰とか出てきて、
「どうやらロジスティック回帰を使うのが定石っぽいけど、プロビットっていう良くわからないのがいっつもくっついて説明されてて困るなぁ」
と思ったりするのは僕だけじゃないはず。そこで自分なりに違いを考えたのでシェアしてみます。
問題１（プロビットが合う）
「ある人の年齢Nを聞いたとき、その人が既婚者か」を確率P(N)で表わすという問題を考えてみます。
結婚という「変化のイベント」について考えると、なんとなく平均結婚年齢あたりにピークがあって、その前後ではなだらかに頻度が少なくなっているイメージがあります。なのでその分布を正規分布だとしましょう。そうすると、年齢Nを聞いたときにその人が結婚してるかどうかは、正規分布の累積分布関数P(N)、すなわちNを与えるとN歳以下で結婚してる確率（N歳現在既婚者である確率）で表せそうです。これがプロビットです。
別の方法として、Pを直接作る代わりに、年齢ごとに既婚者が独身者の何倍いるかを表す関数Q(N)を作ることを考えてみます。もしQが作れれば、既婚者である確率PはQ/(1+Q)ですから簡単です。このQの対数を取ったものlogQを年齢Nと比べることを考えます。このlogQがロジットと呼ばれる値です。対数の底２で考えてみましょう。ある年齢N'のときに既婚者と独身者ちょうど１：１になったとします。このときlogQ=0です。では２：１になるときにlogQはどうなるか？答えは当然１ですね。４：１にならばlogQ=2です。つまり既婚者の比が2倍になる年齢増加がいつも一定というわけです。仮に30歳で既婚者が1:1に、35歳で2:1（全体の6割7分が既婚）ならば、40歳になると4:1つまり8割が既婚者、45歳なら8:1、つまりこれは独身の割合が5歳ごとに半分になっていくというモデルで、なかなか良い気がします。逆に25歳なら1:2つまり既婚者は独身者の半分、20歳なら1:4, 15歳なら1:8（えっ？）、10歳なら（ｒｙとなります。ロジットは正規分布を仮定したプロビットとくらべて年齢が高いほうも低いほうもサチりかたが甘いのが特徴です。
なんとなく違いは分かった気がしますが、印象としてはロジットは単純で分かりやすいもののプロビットのほうが理論としてはイケてるような気がします。理論的背景のあるプロビットのほうがいいんじゃないかと思えてきます。ところが、プロビットはロジットと違って正規分布の累積分布関数の逆関数を計算するのに数値計算アルゴリズムみたいのが必要で、ロジットみたいなlogを取ればOKみたいなお手軽さがないので計算が難しいらしいです。
なのでサチりが甘いという難点がそれほどシビアでない問題の場合はロジットで行こう、みたいなノリらしいです。
けど、ロジットが合う問題というのもあるみたいです。次節。
問題２（ロジットが合う）
「ある人の体重Wを聞いたとき、その人が男かどうか」を確率P(W)であらわすという問題を考えてみます。
まずこれをプロビット的に考えると、これは男性の体重分布と女性の体重分布の境界線をビシッと一つ例えば55キロとかに決めておいて、「各体重における男性の割合」がその５５キロをピークとする正規分布になっていると仮定することになります。これはもう既に仮定自体が矛盾してます。本当は65キロくらいから男性の割合がグッと増えていくのに、女性との境界として設定した５５キロ以降で男性の増分が減っていくと仮定するわけですから。
そういう、それぞれ別の分布を持つ母集団が説明変数をシェアしているだけであって母集団の要素が別の母集団の要素に変化したりしないときにはプロビットは使えなそうです。
これはロジットで考えるほうがマシになりそうです。55キロのときに男性と女性の比が１：１で、50キロだと１：２ならば、60キロにおいては２：１、４０キロなら１：４だし、７０キロなら４：１です。これはなかなか良いような気がします。少なくともプロビットが仮定として崩壊しているのとは雲泥の差です。
まとめ
どうも、説明変数の分布で考えるのがいい気がします。
- 判別したい２値について、説明変数のとりうる値の分布が別々の場合 → ロジット
- 説明変数のとりうる値の分布が共通の場合 → プロビット
- ただしこっちは簡単のためロジットで代用することもまぁアリっちゃアリ
という理解でいいのだろうか。。適材適所で使いわけてるよーって人いましたら教えてくださると嬉しいです。
SciPyで最適化
Deep learningのRBMやAutoEncoderの数式は分かったので、Pythonで実装してみようと思ってSciPyから始めることにします。
非線形最適化と3Dプロットを試してみました。ソースコード。緑の点が初期値、赤の点が最適化結果。緑の４点のぶんだけ最適化関数を試しました。いくつか局所解にハマってますね。
SciPyはニュートン法のような関数の勾配を与える方法だけでなく、勾配情報なしで純粋に関数だけで数値的にやるBFGS法とネルダーミード法というのがあるらしい。BFGSは疑似ニュートン法といって勾配情報を自前で作り上げて動かす方法、ネルダーミード法というのはちょっと変わっていて、三角形をパタパタとひっくり返したり伸び縮みさせたりして尺取虫風に最適な方面に歩いていく方法。
ニューラルネットの重み推定には最尤法を使うそうですが、それはには尤度関数の微分（偏導関数）を求めて勾配法でやるらしいので、SciPyよりも尤度関数の数式自体を記号微分してくれるライブラリtheanoを使うのが便利らしい。
とりあえずもろもろ出来たらまたシェアします。
粒子フィルタを実装してみた
カルマンフィルターを試したので、次のステップとして粒子フィルタをR言語で実装してみました。
ソースはChiral's Gistに置いてあります。
実装
## particle filter implementation by isobe
particle_filter <- function(x0,y,f_noise,f_like,N,M=1) {
tmax <- nrow(y)
D <- length(x0) # == ncol(y)
do_noise <- function(x) {
x1 <- c()
for (i in 1:N) {
for (j in 1:M) {
v <- f_noise(D)
x1<-rbind(x1,x[i,]+v)
}
}
return(x1)
}
do_like <- function(x,t) {
apply(x,1,function(xi) f_like(xi,y[t,]))
}
resampling <- function(x,w) {
# input dim(x) = (N*M,D)
# --> output dim = (D,N)
wsum <- c()
for (i in 1:(N*M)) {
wsum <- c(wsum,sum(w[1:i]))
}
total <- sum(w)
pos <- (1:N) * total/N
r <- runif(1,0,total) # roulette
pos <- (pos+r) %% total
ret <- c()
for (i in 1:N) {
j <- which(wsum>=pos[i])[1]
ret <- cbind(ret,x[j,])
}
return(ret)
}
xx <- list(matrix(x0,D,N))
for (t in 1:tmax) {
x <- xx[[t]] # --> dim(x)=(D,N)
x <- do_noise(t(x)) # --> dim(x)=(N*M,D)
w <- do_like(x,t) # --> length(w)=N*M
xx[[t+1]] <- resampling(x,w) # --> dim(x)=(D,N)
}
return(xx)
}
##### test program #####
test <- function() {
x0 <- c(0,0)
y <- rbind(c(4,4),
c(8,6),
c(6,-1),
c(-2,-5),
c(-8,-9),
c(-6,0),
c(-7,3),
c(-3,6),
c(0,4))
f_like = function(xt,yt) {
return(exp(-sum((yt-xt)**2)))
}
f_noise = function(D) {
rnorm(D,0,3)
}
xx <- particle_filter(x0,y,f_noise=f_noise,f_like=f_like,N=1000)
par(mfrow=c(3,3))
xlim = c(-10,10)
ylim = c(-10,10)
for (t in 1:nrow(y)) {
x <- xx[[t+1]]
plot(y[1:t,1],y[1:t,2],type="b",col="3",xlab="",ylab="",xlim=xlim,ylim=ylim)
par(new=T)
plot(x[1,],x[2,],col="2",xlab=paste("t =",t),ylab="",xlim=xlim,ylim=ylim)
}
}
test()
実行結果
上のコードには２次元平面上でのテストが含まれていて、実行すると以下のようになります。
粒子数N=1000（赤の点）で 時刻10単位分の隠れ状態（緑の点）を推定するという設定です。うまく動いてるようで、良かった。
カルマンフィルタのアドテクへの応用（実践編）
実験は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時間程度でやったのでなかなかやりごたえありました。結果も少し示唆を感じるものが出てきたので面白かったです。
ただ、カルマンフィルタは状態遷移と観測のモデル（行列）のパラメータ数が多く、またモデルの正確さに依拠して解析的にバチンと答えを出す手法なので、ちょっとカタい感じがしました。モデルが正確ならよいですが、パラメータ設定の恣意性は免れません。
それと、今回は対象オーディエンスを集団としましたが、カルマンフィルタで一人一人の興味ベクトルを計算しようとすると頻度の少ない確率事象を線形モデルで表現しずらいという問題があることがわかったので、非線形な粒子フィルタを引き続き検討しようと思いました。
カルマンフィルタの使い方は結構分かったので、他の課題に適用できそうなときがあったらまたやってみたいと思います。
カルマンフィルタのアドテクへの応用（理論編）
カルマンフィルタという数理手法がロケットの姿勢制御等で良く使われています。これをアドテクに応用できそうに思うのでシェアしてみます。
カルマンフィルタと関連手法（知ってる人は飛ばして次へ）
カルマンフィルタは一言でいうと「連続値の隠れ変数についてのモデルが既知でそれが線形で不定性をガウスノイズでくくれるような問題」に適用可能な手法です。
データサイエンティスト的な手法分類としては、線形じゃなくて非線形だったら粒子フィルタというのを使いますが、粒子フィルタはサンプリングといって乱数でシミュレートしてあたりをつけるみたいな方法なのに対して、カルマンフィルタは次の推定値を式一発でビシッと出してくれる（そういうのを「閉形式で解析的に求まっている」みたいな言い方をします）ので線形モデルで近似できそうならば線形モデルを使うことの利点はかなりあります。
隠れ状態が連続値でなく離散値だったら隠れマルコフモデル(HMM)というのを使います。離散値にガウスノイズもくそもないのでHMMの場合は状態遷移のモデル自体を多項分布で表現します*1。
適用方針
カルマンフィルタにおける隠れ状態を「オーディエンス（１名でもいいし、DMPで作ったセグメントとかでもOK）の興味ベクトル」とします。興味というのは例えばアニメが好きだとか、ユニクロというブランドが好きだとか、そういうカテゴリ選好度やブランド選好度だと思ってください。
それとカルマンフィルタは隠れ状態の時間遷移を表す線形変換を指定可能なので「速度」という状態も入れましょう。そうすれば「アニメへの興味が日に日に増している」みたいな各興味要素の変化のスピードを加味できます。
変換行列は対角成分近辺以外はゼロな感じになりそうではありますが対角以外の部分の意味を考えると、競合ブランドの選好度をモデルに組み入れて自社ブランドの速度成分を競合への興味のマイナスとして影響させるようなこともできそうです。
つづいて観測変数ですが「広告への反応ベクトル」とします。おもにクリックやコンバージョン、可視領域の滞在時間とかでもいいかもですね。ベクトルの次元はそれぞれの広告バナーになるので、分析対象とするバナーの種類の数が次元数になります。
そして興味ベクトルから反応ベクトルが導びかれるメカニズムを線形変換（行列）として与えます。この行列は、各広告が反映する興味要素を並べたものになるでしょう。この行列を天下り的に与えないといけないところが難点ですが、とりあえずここは試行錯誤の強い要素だと割り切ることにします。（カルマンフィルタの観測モデル自体をベイズ推定する手法をもし知ってる人いましたら勉強ポイントなど示してもらえるとうれしいです。）
ざっくりとした話になりましたが、以下詳細を見ていきます。
モデルにアドテク変数をあてはめる
まずは興味ベクトルの時間遷移。
は時刻tにおける興味ベクトル、
は上で述べた興味ベクトルの時間遷移です。(時刻tというのは単位時間を表すことにして、例えば1分とか、1時間とか、１日単位とか、分析案件に合わせて設定します。)
の項は上述では触れてませんでしたがロケットでいうとハンドルを切ったりエンジンを噴射したりする制御要素なので、アドテクでいうと「広告を見せた」ということを表すことにして、広告ベクトルと呼びましょう。もし単位時間内に一つも広告を見せなかったらゼロベクトルですし、複数の広告を見せたら複数要素がノンゼロなベクトルです。この広告ベクトルの次元は上述した反応ベクトルの次元と同じです。
は広告ベクトルから興味ベクトルへの変換行列です。とりあえず簡単のため、上で述べた次元が同じである反応ベクトルへの変換行列
(下で出てきます)転置行列（を適当にスケールしたもの）
でいいでしょう。
は多変量ガウスノイズです。これは意味的には一時的な興味の変動で、例えばソーシャルメディアでバズったので一時的にあるテーマへの興味が高まるみたいな外乱や、たまたまモデルの外できっかけが起きて何かに興味を持ったりとかいうのを加味する要素です。ここもいろんな要素を加味しようと思えば上の式のとおりに
を導入して次元数の多い多変量ガウスノイズから射影させることも可能ですが、とりあえず最初は単位行列
として
は興味ベクトルの次元数の分だけ多変量ガウス分布
(
はノイズ生成用の共分散行列)で生成します。
さて、続いて観測モデルです。
は反応ベクトル、
は興味ベクトルから反応ベクトルへの変換行列で、上述した意味合いのものです。
は観測ノイズで
です。
なお、wikipediaを見るとF,H,Q,Rなどの行列はすべて時間に依存する形で書けるみたいですがここではとりあえずこれらは固定としておきます。
予測更新式
前節のようにあてはめた場合の予測式を書き下します。カルマンフィルタが計算してくれるのは、現在時刻の予測値とその誤差の共分散行列
です。
誤差が共分散行列で出てくるというのはイメージ的にいうと、多次元空間上にあるベクトルというのがあって、そのベクトルの矢印のさきっぽを中心として小さなラグビーボールが何らかの方向へ傾いた状態で固定されているような状態です。ラグビーボールの範囲がばらつきの大きさの等高線みたいな感じです。なのでラグビーボールがどんどん大きくなって行ってしまったら困りそうですが、観測変数
で順次補正をかけていくのでうまくいい感じの誤差におさまってラグビーボールが爆発したりしないのでしょう。たぶん。
そしてまず最初に、その初期値と
を与える必要があります。もし初期ベクトルがはっきりわかっていればラグビーボールはゼロ行列でよく、逆に初期ベクトルに不確実性がある場合は適度な大きさのラグビーボールで出発します。
そうしたときの興味ベクトルの予測式が以下になります。
(Tは転置)
ラグビーボールのイメージで考えるとQの項が単純に足し算になっているのがイメージ出来たりして良いですね。
そして反応ベクトルを得たときに興味ベクトルを更新する式が次に来ます。カルマンフィルタ(をはじめとする状態空間モデルにおいてフィルタリングを呼ばれる手法)では観測データの取得前（予測）と取得後（更新）とを分けて考えるのがポイントです。カルマンフィルターはビフォア・アフター、と覚えるといいかもしれません。なので更新ステップとして、上で予測した興味ベクトルと誤差行列を、反応ベクトルを加味した以下の式で上書きします。
まず以下のようにテンポラリ変数を計算しておいて
それを使って
(Iは単位行列)
と更新します。
これで理論が出来上がりました。
続きは実践編で
次回は具体的なデータをでっち上げて実験したものをシェアします。
*1:ちなみに僕が修士課程でやった研究は、統計モデルをPrologで記述するとEM学習やサンプリングなどを自動的にやってくれるPRISMというプログラミング言語を使って将棋倶楽部２４のの棋譜DB約２４万局の表層的なモデリングをするというもので、対局者が強いか弱いかをシステム側で先読みを一切せずに一連の指し手の雰囲気だけから８５％くらいの精度で推定できたのでした。PRISM言語はHMMやベイジアンネットや確率文脈自由文法といった、尤度関数が離散多項分布で表現可能な統計モデルであればPrologプログラムとして表現可能という優れものです。HMMとベイジアンネットを融合したダイナミックベイジアンネットワークとかでもさらっと記述できてEM学習が自動でできます。開発されて１５年以上たつのに全然流行ってないのが残念ですが。
行列計算ライブラリの簡易ベンチマーク
ClojureとF#は行列計算をネイティブコードで実行するライブラリがあるのでそれを使うとどのくらい速いのかをテストしてみました。
3000x3000で要素は[0,1]区間の一様乱数で埋め尽くした行列Mを作ってMとMの逆行列を掛けてトレースを計算する（答えは3000になるはず）という内容です。
条件はMacBookPro late 2013 16GB Memory Win7 Pro x64です。
すべてのライブラリをカバーしてないのと、実行環境によって変化が大きそうなので冒頭での結論提示はしません。参考情報ということで以下でご参照ください。
R (3.0.3 x64)
まずベースラインとしてR言語で。
N <- 3000
pt <- proc.time()
m <- matrix(runif(N*N), nrow=N, ncol=N)
ans <- sum(diag(m %*% solve(m)))
pt <- proc.time()-pt
print(paste("ans =",ans))
print(paste("time = ",pt[3]))何回かやってみて、だいたい平均が43秒でした。ちなみに以下のF#とClojureではもっと速くなる結果が出てますが、F#やClojueの言語そのもので書くとRよりはるかに遅くなります（実際にその言語で書かれたライブラリでもやりました）。
F# (Math.NET Numerics + MKL x64)
open MathNet.Numerics.LinearAlgebra
open MathNet.Numerics.LinearAlgebra.Double
open MathNet.Numerics.Distributions
open System
open System.Diagnostics
MathNet.Numerics.Control.LinearAlgebraProvider <-
new MathNet.Numerics.Algorithms.LinearAlgebra.Mkl.MklLinearAlgebraProvider()
let test n =
let m : Matrix =
upcast DenseMatrix.randomCreate n n (ContinuousUniform (0.0,1.0))
(m * m.Inverse()).Trace()
[<EntryPoint>]
let main argv =
let sw = new Stopwatch()
sw.Start()
let ans = test 3000
sw.Stop()
printfn "%A" ans
printfn "elapsed = %d ms" <| sw.ElapsedMilliseconds
0結果はだいたい10秒程度でした。Rよりだいぶ速い。
.NETは商用のものも多いですがフリーのライブラリはMath.NET Numericsという共通インタフェースにMKLというネイティブライブラリを入れたものがDonSyme著「Expert F# 3.0」によると良いらしいのでそれで試しました。
Clojure (JBLASベースのClatrix)
(ns clojure-matrix-bench.core (:use clatrix.core)) (defn- my-main [] (let [n 3000 m (rand n n)] (trace (* m (i m))))) (defn -main [] (print (time my-main)))
結果はだいたい3〜5秒でした（計測ミス）47〜50秒でした。
Clojureもcore.matrixという共通インタフェースが定義されててそれを各ネイティブラッパーが実装提供する構造になってます。
まとめ
僕のPCはハイパースレッドつき2コアなので物理スレッドは４個ですが、上記のプログラムのいずれも実行中はWindowsのCPUモニタを見るとどうも逆行列計算では1/4しか使っていないっぽかった。（一部、traceのところで５０％に行く時があったので単純な和の場合は並列演算を利かすのでしょう。逆行列計算は途中がモノイドにならないということなのかな）。あとメモリは１億足らずの要素数だからか、全然余裕でした。
というわけでJBlasベースのclatrixがIntel MKLより速かったのがちょっと意外。（計測ミス） Intel MKLは速いですね。
ソースは(こちら)。後でほかのライブラリ（Parallel ColtとEJML）や、ほかの言語のPythonのNumPyとかJulia(たしかBLASのはず)とかでも試してみます。
なにか意味のある結果を出しました的なエントリじゃなくて恐縮ですが、「関数型言語＋データ分析関連ライブラリ」での使える武器を増やしまくりたいという趣旨でもあります。
あと、乱数生成（特にR）が遅い可能性あるという指摘をいただいたので行列演算のみで計測しなおしましたがいずれも1%も変化しませんでした。ただそれは一様乱数だからで、正規乱数とかはパフォーマンスに気を付けたほうがいいかもです。
