スプライン曲線で簡単お手軽ベイズ推定

R Advent Calendar 2013の11日担当記事です。
日付を一日勘違いしてしまい一日遅れになってしまいました。すみません。

テーマは、スプライン曲線で自由に設定した確率分布を使ってベイズ推定をやってみよう、というものです。ガウス分布もベータ分布も出てきません。え、共益分布？なにそれ、尤度関数にマッチする事前分布を設定すればパラメータの更新だけで事後分布出てきますよだって？ノンノン、本記事はそんな小難しいこと知らなくてもベイズ推定のキモは試せるぜ、という話です。

さて能書きはこれくらいにして本題。

R言語でスプライン曲線

飛び飛びの点を用意してsmooth.splineというの関数に与えると、スプライン補間した関数（に相当するオブジェクト）を返してくれます。それに任意のxの点を与えると対応するyの点を返してくれます。
例えば

xy<- rbind(
  c(0,10),
  c(20,100),
  c(50,50),
  c(80,30),
  c(100,10)
)

sp <- smooth.spline(xy[,1],xy[,2],df=5)
d <- predict(sp, seq(1,100,5))

lim <- c(0,100)
plot(xy[,1],xy[,2],xlim=lim,ylim=lim,xlab='x',ylab='y')
par(new=T)
plot(d$x,d$y,pch='',xlim=lim,ylim=lim,xlab='',ylab='')
lines(d$x,d$y)

dfというパラメータは自由度です。点の数と一緒にしておけば大体OK。
グラフはこんな感じになります。

これが確率分布なら分布の形を自由に設定することが出来るというわけです。

この分布からサンプリングする方法は色々あると思いますが、お手軽にやるために離散で考えます。離散多項分布の関数rmultinomというのを使うことにします。上のソースではsmooth.splineで返されるスプライン曲線は任意のxに対するyを計算できるので連続な感じがしていますが、predictはそこから離散的な値をピックアップしていることになります。そのピックアップしたyが多項分布rmultinomのパラメータです。

general_sample <- function(d, size) {
  d$y[which(d$y<0)] <- 0
  r<-rmultinom(1,size,d$y) # d$yの正規化は自動でやってくれる
  rr<-c()
  while (any(r>0)) { 
    i<-which(r>0); 
    rr<-c(rr,i); 
    r[i]<-r[i]-1; 
  }
  return(d$x[rr])
}
hist(general_sample(d,10000))

rmultinomは重みyのインデックスで集計した頻度を行列で返してきますが、今欲しいのはxのサンプリングなのでそうなるよう加工する関数がgeneral_sampleになります。10000個サンプリングしたヒストグラムは以下。

いい感じですね。これで任意の（離散）分布をお手軽に生成する方法が出来ました。

ベイズ推定の問題設定

では例題をでっち上げてモデリングしてみましょう。ある有料課金サービス（ソシャゲーでもビデオオンデマンドでも何でも）の月額課金額を観測変数とし、ユーザの可処分所得を推定するという問題を考えてみます。iさんの課金額をMi, 可処分所得をDiとします。

まず尤度関数（=P(Mi|Di)）を考えます。例えば可処分所得が10万円の人だったら1,2万はまぁまぁのサービスでも使うでしょうけど5,6万は受けるサービスの内容次第です。それが可処分所得5千円の人だったら、千円くらいまでならそこそこのサービスでも使いそうですが、4000円課金するのは厳しそうです。つまり可処分所得のうちの何パーセントかが問題になりそうです。そこで可処分所得のうちのパーセンテージごとの課金傾向分布は固定にします。とりあえず可処分所得の２０％をピークとする分布を仮定してみます。

実は、上で描いたグラフがその分布を表していました。

かなり恣意的な感じがしますが、ベイズ推定なんてそんなもんです。事前分布も仮説、尤度関数も仮説、仮説仮説でPDCA回して本当に効果ありましたか？というメソッドです。という意味でスプライン曲線を使うとこのあたり自由に色々出来てよいです。

続いて事前分布（=P(Di)）ですが、とりあえず可処分所得は1000円から10万円まで一様に分布していることにします。これも分布の作り方は上記と同様です。他の設定を試したければ後でスプラインの点を色々変えればいいでしょう。

そして事後分布（=P(Di|Mi)）がどうなるかですが、上のソースにもあるように正規化はrmultinom関数が自動的にやってくれるので事前分布と尤度関数を単純に掛け算すればOKです。だた観測データである課金額は複数与えて計算したいので、各消費者の課金は独立でしょうから「P(Di|Mi)∝P(Di)*P(Mi|Di)」をすべてのiについて掛け算すればOKです。というわけでRの関数は以下のようになります。

like_base <- rbind(
  c(0,10), c(20,100), c(50,50), c(80,30), c(100,10)
)

like_sp <- function(di) {
  spline(cbind(di*like_base[,1]/100,like_base[,2]),4)
}

pre_sp <- spline(rbind(
  c(1000,10),
  c(5000,10),
  c(10000,10),
  c(100000,10)
),4)

pre_di <- predict(pre_sp,seq(1000,100000,1000))

post_di <- function(paid) {
  d <- pre_di
  for (i in 1:length(d$x)) {
    x <- d$x[i]
    y <- d$y[i]
    w <- y*predict(like_sp(x),paid)$y
    w[which(w<1)] <- 1
    d$y[i] <- exp(sum(log(w))) # multiply for all w
  }
  return(d)
}

尤度関数like_spは上述した可処分所得のパーセンテージごとの分布を元に、引数で与えられた可処分所得に応じて横方向に伸び縮みさせた曲線を返しています。事後分布post_diは観測データである課金額paidを受け取って各可処分所得に応じたpaidの出現度と事前の重みを掛けます。（尤度関数はパラメータと観測データの2変数関数なので、イメージとしては3次元グラフを観測データごとにスライスしたものを掛け合わせている感じになります。）

テスト

試しにいくつか観測データを与えて事後分布を描いてみましょう。1000〜3000円くらいが4名と、1万円ブッコむ大人が1名という状況をインプット。

d <- post_di(c(3000,2000,1000,10000,2500))
plot(d$x,d$y)
lines(d$x,d$y)

出ました。ちょっと指数表示になって見づらい（こういうデフォルトの動作が？？？なのが僕のRの一番嫌いなところ）ですが、だいたい可処分所得1.7万くらいにピークがありますね。形も尤度関数で設定したパーセンテージのグラフに似てます。ただ仮説が恣意的とはいえ、5名のデータを与えた結果がその恣意性から定量的に計算されるのはやっぱりベイズの強力なところですね。

ちょっと設定を変えてやってみます。同じ尤度関数、観測データもそのままで、事前分布をもう少しそれっぽい（サラリーマンのこずかい平均3万円を意識して）1〜3万円くらいに山があって以後10万までなだらかに減っていくような分布を仮定してみます。

これで事後分布を描くと以下のようになりました。

最初にやったフラットな事前分布の時よりも可処分所得のピークが2万円くらいと少し高くなっていて、ピークの山もなだらかになっています。事前分布がフラットだとピークがとがるのは事前分布が一様な場合に事後分布は最尤推定っぽくなるという話と一緒だと思います。なだらかになったのは、事前分布がフラットなときと比べて最尤なポイントの近辺で事前分布が手厚いのでその名残りがベイズ更新後に現れているということでしょうかね。あと、フラットな事前分布の時と比べてすそ野が削られているのも事前分布の形の影響ですね。

ここまで重みのグラフをプロットする形で話してきましたが、もちろん最初に述べたサンプリング関数を使ってこのモデルからサンプリングも出来ます。

あとオーバーフローの問題があるので注意は必要です。正規化をせずに重みを掛けていくので観測データが多いときはいくつかに分けて、一定数の観測データを食わせた事後分布をいったんいくつかの代表点によるスプライン関数に落とす、というのを繰り返すことになると思います。

ソース全体はGistにおいてありますので興味ありましたらどうぞ。

まとめ

スプライン曲線を使って任意の分布をお手軽に作ってベイズ推定するという話をしました。ただこの方法は多次元でパラメータ数が多くなってくると難しくなります。尤度関数の計算が確率変数の空間全体から観測データの軸でスライスして重みを足さないといけないからです。
ですが上述の例題のように1次元であれば複雑な仮定を置いて定量的に検証したりは出来ますし、多次元でも値の少ない離散値（年齢や性別とか）なら対応できそうです。上述の例題では可処分所得から何ペーセントを課金するかを固定してますがそれを変えたりとか、1次元でも出来ることは多いと思います。

ですが、まぁこれはお遊びですかね(^^)