数値計算的にやろう あっちかこっちか

• ベータ分布を用いた理論的計算を前の記事で書いた
• 既知分布を使えないこともある
• そんなときに、うまく数値計算的に同じことができると便利

# 数値計算的にやろう
# 既存の関数などが使えないが、２次元空間に格子状の確率値が得られたとして、同じことを実行する
# qq.both.pre, qq.sep.pre, qq.both, qq.sepが２次元空間格子点値の行列として与えられているとする
# それぞれの値が、微小単位正方形における確率を表すようにするためには
# 全値が代表している単位正方形のサイズが等しいので、総和が１になるように調整すればよい
qq.both.pre.a <- qq.both.pre/sum(qq.both.pre)
qq.sep.pre.a <- qq.sep.pre/sum(qq.sep.pre)
qq.both.a <- qq.both/sum(qq.both)
qq.sep.a <- qq.sep/sum(qq.sep)
# これを用いて、事前、事後の仮説別尤度の変化比を求める
# 男女同じの仮説のそれは、男女別の行列を用いて求める
pr.both.pre <- sum(diag(qq.sep.pre.a)) # 男女別の行列の対角成分を用いている
pr.sep.pre <- sum(qq.sep.pre.a)
pr.both.post <- sum(diag(qq.sep.a)) # 男女別の行列の対角成分を用いている
pr.sep.post <- sum(qq.sep.a)
# 事前事後での、「総尤度」の変化から、両者の尤度変化の比を出す
r.both <- pr.both.post/pr.both.pre
r.sep <- pr.sep.post/pr.sep.pre
# pre.pからpost.pを出す
post.p.a <- pre.p*c(r.sep,r.both)
post.p.a <- post.p.a/sum(post.p.a)
# 数値計算的な結果と前出の『理論値』を比べる
post.p.a
post.p

# これを用いて「女に想定する生起確率」を計算しよう
p.women.pre.a <- apply(qq.sep.pre.a,2,sum)
p.women.post.a <- apply(qq.sep.a,2,sum)
p.both.pre.a <- apply(qq.both.pre.a,2,sum)
p.both.post.a <- apply(qq.both.a,2,sum)

est.women.pre.a <- p.women.pre.a * pre.p[1] + p.both.pre.a * pre.p[2]
est.women.post.a <- p.women.post.a * post.p[1] + p.both.post.a * post.p[2]

# 格子サイズに関する補正
# 先に、総和が１になるように補正し、その結果として出てきている値であるから
# セルの数(この場合は１次元空間なので、セルの数はlength(p)で再補正し、
# 理論値と合うようにする
est.women.pre.a <- est.women.pre.a *length(p)
est.women.post.a <- est.women.post.a *length(p)
# 格子サイズに関する補正
# 先に、総和が１になるように補正し、その結果として出てきている値であるから
# セルの数(この場合は１次元空間なので、セルの数はlength(p)で再補正し、
# 理論値と合うようにする
est.women.pre.a <- est.women.pre.a *length(p)
est.women.post.a <- est.women.post.a *length(p)
# 理論値と数値計算経由の値とを比較する
par(mfcol=c(1,2))
plot(est.women,est.women.post.a)
abline(0,1,col=2)
matplot(p,cbind(est.women,est.women.post.a),type="l",main="事後")
legend(max(p)*0.6,max(p.men)*0.9,c("est-women","数値計算"),lty=1:2,col=1:2)
par(mfcol=c(1,1))