- パッケージ 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ではランダムにグラフを作ってくれる。上で、5要素のデータを扱ったので、5要素のグラフを100個作らせてみる。作成にオプションがあるのは、いつものことだが、デフォルトで作るとすれば以下の通り
random.graph(LETTERS[1:5], num = 100)
-
- グラフの形の推定のためには、「グラフのよしあし」の基準が必要。これにも基準があるが・・・。推定のグラフ(res)を、データ(d)にて、採点する
- ただし、場合によっては、推定グラフの辺の向きが決まっていないものがある。その場合は採点ができないので、それをどちらかにする必要がある(そのための関数が用意されているし、入力することもできる。
scoregs<-score(res, d)
-
- ランダムに作ったグラフも、データ(d)にて採点する。ランダムグラフのうち、4個だけ、描いておく。
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)