確率的グラフィカルモデル - 驚異のアニヲタ社会復帰の予備

読んだ。

作者: 鈴木譲,植野真臣,黒木学,清水昌平,湊真一,石畠正和,樺島祥介,田中和之,本村陽一,玉田嘉紀
出版社/メーカー: 共立出版
発売日: 2016/07/23
メディア: 単行本
この商品を含むブログ (2件) を見る

COI：ラボにあった。編集者とは名刺交換したことがある程度。
　
ベイジアンネットワークのようなグラフィカルネットワークについての話はあまり詳しくないので読んだ。複数人が各章を好き勝手に書いているが、用語や記法は統一されているのでそこまで混乱しない。
統計力学とグラフィカルモデルはほとんど知らない分野だったので意味不明なところがほとんどだったが、それ以外はグラフィカルモデルの基礎と、ゲノム解析に応用した例もあり、なかなかよかった。
　
用語や定義を覚えていないので写経しておく。
d分離という概念が重要である。A→B→C という逐次結合の時、B の状態が判明してしまったら、AとCは独立になってしまう。このとき、AとCは「Bを所与としてd分離」という。B のように、あるひとつの状態がわかってしまうことを「その変数がインスタンス化された」という。
　
A→{B,C,...}のときは分岐結合、{B,C,...}→A のときは合流結合というが、合流結合のときは、変数{B,C,...} はAもしくはAの子孫について証拠を得ていないときに「Aを介してd分離である」というらしい。
まとめると、因果モデルにおける2つの変数A, Bは、AとBのすべてのパスに存在する以下のような変数V があるとき、d分離という。

逐次結合もしくは分岐結合で、Vがインスタンス化されているとき

または

合流結合で、VもしくはVの子孫がインスタンス化されていないとき

AとBがd分離でないとき、d結合と呼ばれる。
　
本ではインスタンス化された因果モデルの例が出ている。エビデンスを得てインスタンス化されたノードを灰色にしている。
まず、エビデンスが何もない場合は、ノードA の情報はノードB には伝搬しないため、d分離である。
エビデンス1 のみが得られた場合は、Aの情報はBの子ノードまでは伝搬されるが、B には伝搬されないのでd分離されている。
エビデンス1と2を同時に得たときは、A の情報は合流結合のルールを2回繰り返してBに伝搬するため、d結合である。
エビデンス1, 2, 3 が同時に得られると、エビデンス3 がエビデンス2 の合流結合によるd分離を阻止するため、Bに情報が伝搬せず、d分離となる。
らしい。

library(igraph)
el <- cbind(c(1,2,3,3,4,5,6,6,7,8),c(5,5,6,7,7,8,8,9,9,10))
g <- graph_from_edgelist(el)

V(g)$label <- replace(rep(NA, length(as_ids(V(g)))), c(1,4,6,9,10), c("A", "B", 3:1))
V(g)$label.cex <- 2
V(g)$label.font <- 2
V(g)$color <- replace(rep("white", length(as_ids(V(g)))), c(6,9,10), "grey")
E(g)$width <- 2
E(g)$arrow.width <- 2
E(g)$color <- "black"

lay0 <- do.call(rbind, mapply(function(z) cbind(seq((4:1)[z])+0.5*(z-1), -sqrt(2)*(z-1)), 1:4))
par(mar=rep(0, 4))
plot(g, layout=lay0)

　
バックドア基準とは、DAG であるグラフG において、XはY の非子孫であるとする。このとき、

X からZ へ任意の要素の有向道がない
G よりX から出る矢印をすべて除いたグラフにおいて、Z が X と Y をd分離する

とき、変数集合Z は(X, Y) についてバックドア基準を満たす、という。
　
いま、図のような因果モデルを考える。X3→X4 の因果関係について、どんな変数がバックドア基準を満たすかどうかやってみよう。

　
Z がX2 である場合を想定すると、Z=X2, (X, Y)=(X3, X4) である。
X3 からX2 への有向道はない(1)、G よりX3 から出るすべての矢印を取り除いた下の図で、X2 は(X3, X4)の分岐結合であり、X2 はX3 とX4 をd分離している。
これは{X1, X2}の集合についても言える。これにより、X2 は(X3, X4)についてバックドア基準を満たしており、X2 を観測することでX3 からX4 への因果的効果は識別可能になる。

　
フロントドア基準とは、DAG であるグラフG において、XはY の非子孫であるとする。このとき、

X からY への任意の有向道上にZ の要素が存在する。
G よりX から出る矢印をすべて除いたグラフにおいて、空集合はX とZ の任意の要素をd分離する。
G よりZ の任意の要素から出る矢印をすべて除いたグラフにおいて、X はZ の任意の要素とY をd分離する

下ふたつの要件はそれぞれ

G において、空集合は(X, Z) についてバックドア基準を満たす。
G において、X は(Z, Y) についてバックドア基準を満たす。

のと同じことらしい。
上のグラフからX2→X4 の有向道がないグラフについて、X2→X4 の因果的効果を定式化してみる。つまり、Z=X3, (X, Y)=(X2, X4) である。
X2 からX4 にいたる道にX3 があるのは見ての通りで、X2 から矢印を除去するとX2(=X) とX3(=Z) は空集合でd分離される。最後に、X3(=Z) から矢印を除去すると、X2(=X) はX3 とX4(=Y) をd分離していることがわかる。よって、X3 は(X2, X4) についてフロントドア基準を満たす、といい、X3 を観測することでX2→X4 への因果的効果は識別可能となる。

el <- matrix(paste0("X", c(1,2,1,4,2,4,2,3,3,4)), nc=2, byrow=TRUE)
g0 <- graph_from_edgelist(el)
lay0 <- cbind(c(267,109,387,182), c(200,117, 0, 2))
V(g0)$size <- 40
V(g0)$label.cex <- 3
V(g0)$label.font <- 2
E(g0)$color <- "black"
E(g0)$width <- 5
E(g0)$arrow.width <- 2
plot(g0, layout=lay0)

E(g0)$color[grep("^X3", as_ids(E(g0)))] <- "grey"
E(g0)$color[grep("^X3|X2\\|X4", as_ids(E(g0)))] <- "grey"

　
条件付き操作変数法とは、因果グラフG において、X はY の非子孫とする。このとき、

Z とW に含まれる任意の頂点はX およびY の子孫でない。
G よりX から出る矢印を除いたグラフにおいて、W はZ とY をd分離するが、Z とX をd分離しない。

を満たす変数Z は、変数集合W を与えたときの(X, Y) に対する条件付き操作変数法という。
　
構造方程式と非ガウス性
変数がX1, X2 のふたつあるとき、因果の方向はX1→X2 かX2→X1 である。データ行列が与えられたときにこの因果関係を知りたいのだが、十分条件として次のふたつが成り立てばよい。

e1 あるいはe2 が非ガウス分布に従う。
e1 とe2 が独立である。

e1 とe2 はそれぞれX1, X2 を生成する連続な誤差変数で、図示するとこうなる。いま、係数 $\beta$ をとって2つのモデル
$\begin{align}x_1&=&e_1\\x_2&=\beta_{21}x_1+e_2\end{align}$ もしくは $\begin{align}x_1&=&\beta_{12}x_2+e_1\\x_2&=e_2\end{align}$ のどちらかを判定したい。

　
非ガウス分布であればいい、という緩い仮定でよく、セミパラメトリックと言われる。
rで標本数 $n$ だけを設定したらランダムサンプリングをしてくれる分布を選んでやってみると、確かに正規分布(norm)以外の分布では、左右を見比べると分布の形に差があるので、モデルの判別ができる。
(ロジスティック分布は $n$ が小さいと判別こんなんだが、増やすと差が歴然としてくるようになった。)

dens <- c("p", "r", "q", "d")
fs <- mapply(function(z) gsub(paste0("^", z), "", apropos(paste0("^", z))), dens)
fs <- intersect(intersect(intersect(fs$q, fs$p), fs$r), fs$d)
fs <- fs[c(13,3,5,10,11,17)]

n <- 20000
hoge <- mapply(function(z) matrix(eval(parse(text=paste0("r", z, "(", n*2, ")"))), nc=2), fs, SIMPLIFY=FALSE) # 分布からランダムサンプリング

b <- 0.8 # 線型結合の係数
m1 <- m2 <- array(0, c(n, 2, length(fs)))
for(i in seq(fs)){
  m1[,,i] <- cbind(hoge[[i]][,1], b*hoge[[i]][,1] + hoge[[i]][,2])
  m2[,,i] <- cbind(b*hoge[[i]][,2] + hoge[[i]][,1], hoge[[i]][,2])
}

par(mfrow=c(length(fs), 2), cex.main=2, pch=16)
for(i in seq(fs)){
  plot(m1[,,i], xlab="X1", ylab="X2", main=paste0(fs[i], ": ", "X1 ", "to ", "X2 ", "model"))
  plot(m2[,,i], xlab="X1", ylab="X2", main=paste0(fs[i], ": ", "X2 ", "to ", "X1 ", "model"))
}