2013-09-06

あっちかこっちか

決断診断ベイジアンネットワークベータ分布ロジスティック回帰射影多様体微分幾何カルマンフィルタ

これは下書き
昨日の話題の脇道
要因Xのありなしとともに、ある事象Yの生起の有無の集計をとったら

	Y(+)	Y(-)
X(+)	a	b
X(-)	c	d

だったという
さて。新しいサンプルがX(+)だったときに、Yの生起確率に何を思うのか…
X(+),X(-)で生起確率が違うのなら、(a,b)を用いてベータ分布
X(+),X(-)で生起確率が同じなら、(a+c,b+d)を用いてベータ分布
問題は、X(+),X(-)で生起確率が同じなのか違うのかについて判然としないこと
判然とはしないけれど、事前確率として、X(+),X(-)で生起確率が同じ・違うが0.5,0.5だという「気持ち」はあるとする、『ベイズ』でやっているから
X(+),X(-)で違うときは
- $(p_+,p_-)$ という要因別の生起確率の組が作る2次元正方形空間の座標の尤度(!)を $B(a+1,b+1)(p_+) B(c+1,d+1)(p_-)$ としておくと
- $\int_0^1 B(a+1,b+1)(p_+) B(c+1,d+1)(p_-) dp_+ dp_- =1$ となる
- X(+)の生起確率だけが気になれば、ただの $B(a+1,b+1)$ になるのは、2次元正方形空間をX(-)について積分して $\int_0^1 B(a+1,b+1)(p_+) B(c+1,d+1)(p_-) dp_- = B(a+1,b+1)(p_+)$ としてもよいし、 $B(a+1,b+1)(p_+)$ だよ考えなくても、ということにしてもよい
X(+),X(-)で同じときには
- $p_+ = p_- = p$ なので
- $k B(a+1,b+1)(p) B(c+1,d+1)(p)$ として
- $\int_0^1 B(a+1,b+1)(p) B(c+1,d+1)(p) dp =K$ になるような $k,K$ を考えてやればよい
- ここで $k,K$ と２つの変数がある $K$ はX(+),X(-)で生起確率が同じであるという仮説の「事後確率」に当たる
- 「空間」で考えれば、2次元正方形空間の対角線部分に相当する1次元線分空間の積分
２つの確率密度関数は方や１次元、方や２次元なので、それが積分してある値(スカラー)になる、とは言ってもちょっと面倒くさい
『どのくらいの割合でおきるか』と考えると土俵が同じになるので、そうしてみよう
$(p_+,p_-),(p_+ + \delta,p_- + \delta)$ の2点を対角点とする微小正方形の生起確率と
$(p),(p+\delta)$ の2点を両端とする線分〜 $(p_+=p,p_-=p),(p+\delta,p+\delta)$ の2点を対角点とする微小正方形の生起確率との多寡を比べることにしよう
前者は $B(a+1,b+1)(p_+) B(c+1,d+1)(p_-)\times \delta^2$
後者は $B(a+1,b+1)(p_+=p) B(c+1,d+1)(p_-=p)\times \delta^2$
これをもとにの確率をX(+),X(-)が違う・同じの２仮説別に計算すると
- $\int_0^1 \int_P^{P+\delta} B(a+1,b+1)(p_+=P) B(c+1,d+1)(p_-) dp_+ d p_-$
- $\int_P^{P+\delta} \int_P^{P+\delta} B(a+1,b+1)(p_+=p) B(c+1,d+1)(p_-=p) d p_+ d p_-$
これを数値計算風にするべく、離散的に計算するなら
- $\sum_{i=1}^n B(a+1,b+1)(p_+=P) B(c+1,d+1)(p_-=p_i) \times \delta^2$
- $B(a+1,b+1)(p_+=P) B(c+1,d+1)(p_-=P) \times \delta^2$
事前確率のことを考えると、X(+),X(-)が異なるという仮説では、すべての微小正方形 $\delta^2$ が等確率であると考えていたなら、その確率は $\frac{1}{\delta^2}$ 。同様にX(+),X(-)で等しいという仮説では、すべての微小正方形が $\frac{1}{\delta}$ となる。
となって、結局、X(+)の際にそうていするべき、事後の生起確率はこの事前の生起確率を加味して
- $\sum_{i=1}^n B(a+1,b+1)(p_+=P) B(c+1,d+1)(p_-=p_i) \times \delta^2\times \frac{1}{\delta^2}$ + $B(a+1,b+1)(p_+=P) B(c+1,d+1)(p_-=P) \times \delta^2\times \frac{1}{\delta}$
なお、X(+),X(-)に異なる生起確率を想定するのがよいか、想定しないのがよいか、というのを仮説検定するのなら、それぞれのモデルでの最尤推定尤度を比較して、自由度の差(ここでは1)で検定する、ということになるが、今は、検定にはまったく興味がないので、それは放っておく
上のやり方と別の方法でもできそうだが、その方法と上の方法で値が合わない…どっちかが間違っている…

men <- c(1,5)*5
women <- c(4,2)
both <- men + women

p <- seq(from=0,to=1,length=100)

p.men <- dbeta(p,men[1]+1,men[2]+1)
p.women <- dbeta(p,women[1]+1,women[2]+1)
p.both <- dbeta(p,both[1]+1,both[2]+1)

par(mfcol=c(2,2))
#matplot(cbind(p.men,p.women,p.both),type="l")

q.men <- exp(lgamma(sum(men)+1)-sum(lgamma(men+1))+men[1]*log(p)+men[2]*log(1-p))
q.women <- exp(lgamma(sum(women)+1)-sum(lgamma(women+1))+women[1]*log(p)+women[2]*log(1-p))

q.both <- q.men * q.women
qq.sep <- outer(q.men,q.women,"*")
qq.both <- diag(q.both)

#matplot(cbind(c(qq.sep),c(qq.both)),type="l")

pre.p <- c(0.5,0.5)
#qq.sep.st <- qq.sep/sum(qq.sep)*pre.p[1]
#qq.both.st <- qq.both/sum(qq.both)*pre.p[2]
qq.sep.st <- qq.sep/((length(p)-1)^2)
qq.both.st <- qq.both/(length(p)-1)

image(qq.sep.st+qq.both.st)
contour(qq.sep.st+qq.both.st,add=TRUE)


a <- apply(qq.sep.st,2,sum)
b <- apply(qq.both.st,2,sum)
ab <- a+b

q.women.2 <- ab/sum(ab)
p.men.2 <- p.men/sum(p.men)
p.women.2 <- p.women/sum(p.women)
p.both.2 <- p.both/sum(p.both)
matplot(p,cbind(p.men.2,p.women.2,p.both.2,q.women.2),type="l")
legend(max(p)*0.6,max(p.men.2,p.women.2,p.both.2,q.women.2)*0.8,c("men","women","both","combination"),lty=1:4,col=1:4)

# 両仮説からの分布を事前確率重みづけで合算すればよいのかというとそうではない
q.women.3 <- p.women*pre.p[1]+p.both*pre.p[2]
#plot(q.women.2,q.women.3)

matplot(p,cbind(q.women.2/sum(q.women.2),q.women.3/sum(q.women.3)),type="l")
legend(max(p)*0.6,max(q.women.2/sum(q.women.2),q.women.3/sum(q.women.3))*0.9,c("combination","linear sum"),lty=1:2,col=1:2)

par(mfcol=c(1,1))