メモ

N <- 100
x <- seq(from=0,to=1,length=N)
x <- x[c(-1,-N)]
y <- x
X <- c(40,40)
Y <- c(40,5000)
beta.x <- dbeta(x,X[1]+1,X[2]+1,log=TRUE)
beta.y <- dbeta(y,Y[1]+1,Y[2]+1,log=TRUE)
#prev.data <- c(560,1000)
prev.data <- X+Y
p <- x
beta.p <- dbeta(p,prev.data[1]+1,prev.data[2]+1,log=TRUE)

xyp <- expand.grid(x,y,p)

frac.x <- xyp[,1]*xyp[,3]/(xyp[,1]*xyp[,3]+xyp[,2]*(1-xyp[,3]))
tmp <- apply(expand.grid(beta.x,beta.y,beta.p),1,sum)
prob.xy <- exp(tmp-max(tmp))

#plot(frac.x,prob.xy)

plot(density(frac.x,weights=prob.xy/sum(prob.xy)))
points(p,exp(beta.p),type="l",col=2)


abline(v=prev.data[1]/sum(prev.data))
abline(v=(X[1]/sum(X)*prev.data[1]/sum(prev.data))/(X[1]/sum(X)*prev.data[1]/sum(prev.data)+Y[1]/sum(Y)*(1-prev.data[1]/sum(prev.data))),col=2)