2010-03-19

bnlearnでベイジアンネットワーク

パッケージ　bnlearn のホームページとCRAN
ベイジアンネットワークの概要
- データがある
- 要素をAcyclic directed graph状にして、一番よい形を見つけ、その形の辺に適当な値を見出す
bnlearnでなぞる
- データはdata.frame。こんな風にして、関係のあるデータを作ってやってもよい。

Ns<-10000

A<-rpois(Ns,10)
B<-rpois(Ns,10)
C<-rpois(Ns,10)
D<-rpois(Ns,10)
E<-rpois(Ns,10)

B[1:(4*Ns/5)]<-A[1:(4*Ns/5)]
C[1:(3*Ns/5)]<-A[1:(3*Ns/5)]
D[1:(2*Ns/5)]<-A[1:(2*Ns/5)]
E[1:(1*Ns/5)]<-A[1:(1*Ns/5)]

d<-as.data.frame(cbind(A,B,C,D,E))
# 実数と認識させる
d <- d+0.0

- グラフの形の推定
  - グラフの形の推定は、「しらみつぶし」にはそぐわない
  - しらみつぶしをせずに探す方法はネットワーク探索の大事な部分。いろいろ方法がある
  - bnlearnもいくつかの方法を実装している。その内容に触れることは、この記事の目的ではない。
    - costraint-based algorithms:
      - Grow-Shrink (GS)
      - Incremental Association Markov Blanket (IAMB)
      - Fast Incremental Association (Fast-IAMB)
      - Interleaved Incremental Association (Inter-IAMB)
    - the following score-based algorithms:
      - Hill Climbing (HC)
      - Tabu Search (Tabu)
    - and the following hybrid algorithms:
      - Max-Min Hill Climbing (MMHC)
      - General 2-Phase Restricted Maximization (RSMAX2)
  - 方法のひとつgsで推定して、グラフを描く

library(bnlearn)
res<-gs(d)
plot(res)

- グラフの探索アルゴリズムがあるのはよいとして、それがなければ、「全部試すか、ランダムに試す」というのが、よくやること
  - bnlearnではランダムにグラフを作ってくれる。上で、５要素のデータを扱ったので、５要素のグラフを100個作らせてみる。作成にオプションがあるのは、いつものことだが、デフォルトで作るとすれば以下の通り

random.graph(LETTERS[1:5], num = 100)

- グラフの形の推定のためには、「グラフのよしあし」の基準が必要。これにも基準があるが・・・。推定のグラフ(res)を、データ(d)にて、採点する
  - ただし、場合によっては、推定グラフの辺の向きが決まっていないものがある。その場合は採点ができないので、それをどちらかにする必要がある(そのための関数が用意されているし、入力することもできる。

scoregs<-score(res, d)

- ランダムに作ったグラフも、データ(d)にて採点する。ランダムグラフのうち、４個だけ、描いておく。

Niter<-100
Ndraw<-4
par(mfcol=c(2,2))
scores<-rep(0,Niter)
for(i in 1:Niter){
 randomGraph<-random.graph(LETTERS[1:5], num = 1)
 if(i<=Ndraw){
  plot(randomGraph)
 }

 scores[i]<-score(randomGraph,d)
 #randomGraphs<-c(randomGraphs,randomGraph)
}

- 採点結果をプロットして、推定されたグラフの採点結果(水平線)と較べる(掲載図)
  - 十分なサンプル数で、まともな「因果関係」が潜んでいれば、推定グラフが最高点を出す。

par(mfcol=c(1,1))
ylim<-c(min(scores,scoregs),max(scores,scoregs))
plot(scores,ylim=ylim)
abline(h=scoregs)

- 要素間の関係に数値を与える(辺に値をつける)
  - 離散的な場合は、条件付確率
  - 連続的な場合には、回帰係数

bn.fit(res,d)

例題２

Ns<-10000

A<-rpois(Ns,10)
B<-rpois(Ns,10)
C<-A+B
D<-A+B+C
E<-A^2
F<-rpois(Ns,30)
G<-rnorm(Ns)
H<-F+G
d<-as.data.frame(cbind(A,B,C,D,E,F,G,H))
d <- d + rnorm(length(d))*2
res<-gs(d)
plot(res)