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 \times B(a+1,b+1)(p) B(c+1,d+1)(p)$ として
- $\int_0^1 k \times B(a+1,b+1)(p) B(c+1,d+1)(p) dp =K$ になるような $k,K$ を考えてやればよい
- 実際、 $B(a+1,b+1)(p) B(c+1,d+1)(p) \propto B(a+c+1,b+d+1)(p)$ であり、 $\int_0^1 B(i,j)(p) dp =1$ であるので、
- $k \times B(a+1,b+1)(p) B(c+1,d+1)(p) = K \times B(a+c+1,b+d+1)(p)$
- $K = \frac{B(a+1,b+1)(p) B(c+1,d+1)(p)}{B(a+c+1,b+d+1)(p)}$ が任意の $p$ について等しい値であって、それはベータ分布の定義式から
- $K=\frac{(a+b+c+d+1)!}{(a+c)!(b+d)!}\frac{a!b!c!d!}{(a+b+1)!(c+d+1)!}$ で求められる

# 男女別に集計
men <- c(20,4)
women <- c(2,4)
# 男女合算集計
both <- men + women
# 計算対象とする生起確率(0,1を除いて、(0,1)で均等配置)
p <- seq(from=0,to=1,length=50)
p <- p[-length(p)]
p <- p[-1]
# それぞれ、ベータ分布
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)
# ベータ分布の「由来」はこの尤度
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))
# 男女合算という仮説で出したベータ分布p.bothは「本当の尤度」ではなく、男女で生起確率が等しいという『狭い仮説空間』の中で生起確率の総和が１になるように調整された値であって
# 男女別もありうる、という仮説では、その確率は下がる
# その「実際に下がった値」を計算しておく
q.both <- p.men * p.women
# 男女別もありを想定した仮説空間にあって、男女同じ仮説に限定したら、２次元正方形仮説空間のうち、対角線に相当する１次元線分にのみ正の値のある確率密度分布(様)のものがある(２次元積分が通常の意味でうまく行かないので『確率密度分布(様)』と呼んでおく
qq.both <- diag(q.both)
# 男女別を前提とした確率密度分布は、単純に、男の生起確率と女の生起確率の『積な分布』
qq.sep <- outer(p.men,p.women,"*")

# 男女同じと限定したときのベータ分布と、男女別もありを仮定した上で、「たまたま男女が同じ生起確率」というときの、その仮説位置(対角線)の確率密度分布の値は、定数倍の関係
my.k <- (q.both/p.both)
my.k
# この定数は(a,b)という観察のときのベータ分布が(a+b+1)!/(a!b!)p^a(1-p)^bであることから次のようにしても出る
my.k.2 <- exp(-(lgamma(sum(both)+2)+sum(lgamma(c(men,women)+1)))+sum(lgamma(c(sum(men),sum(women))+2)+lgamma(both+1)))
my.k.2

# 男女が等しいという仮説と男女で異なる(異なってもよい)という仮説とについて、事前確率を持っていたとする。その事前確率
pre.p <- c(0.5,0.5)
# 事前確率(データが無かったときの確率)は
p.men.pre <- dbeta(p,1,1)
p.women.pre <- dbeta(p,1,1)
p.both.pre <- dbeta(p,1,1)
# であって、これを２次元仮説空間に展開しても
qq.both.pre <- diag(p.men.pre*p.women.pre)
qq.sep.pre <- outer(p.men.pre,p.women.pre,"*")

# データを受け取った後で、２つの仮説の尤度がどう変化したかを考える
# 男女別という仮説が、qq.sep.pre -> qq.sep に変化したと考える
# 両者の空間全体の積分はともに１である(値に変化はないが、『相対的には変化している』と考える
# 他方、男女同じという仮説では、qq.both.pre -> qq.both に変化したと考える
# qq.both.preの積分(１次元線分上の積分)は１
# qq.bothの積分(１次元線分上の積分)は、p.bothというベータ分布(これは積分が１）と形が同じで、my.k定数倍である分布の積分はmy.kになっていると考える
# 男女別仮説は1->1、男女同仮説は1->my.k
# これを利用して、男女が異なる・同じの２仮説の相対尤度は
# pre.p[1] * 1/1,pre.p[2]*my.k/1へと変化したと考えて
post.p <- pre.p * c(1,my.k.2)
post.p <- post.p/sum(post.p)
# これを利用すれば、女に想定する生起確率は
est.women <- p.women * post.p[1] + p.both * post.p[2]
# 同じやり方での事前確率は
est.women.pre <- p.women.pre * pre.p[1] + p.both.pre * pre.p[2]
par(mfcol=c(2,2))
matplot(p,cbind(est.women.pre,est.women),type="l",main="女予測、事前事後")
# 男女別、男女同じをpre postで分けて示せば
matplot(p,cbind(p.men.pre,p.women.pre,p.both.pre,est.women.pre),type="l",main="事前")
matplot(p,cbind(p.men,p.women,p.both,est.women),type="l",main="事後")
legend(max(p)*0.6,max(p.men)*0.9,c("men","women","both","est-women"),lty=1:4,col=1:4)

par(mfcol=c(1,1))