p<-0.5
n<-c(100,1000,10000,50000)
n*p*(1-p)

library(wesanderson)
mycolors<-wes_palette("GrandBudapest1",4)

## prior is 0.1, 0.1
par(mfrow=c(1,3), pty='s')
prange<-seq(0.4, 0.6, 0.001)
n<-c(100,1000,10000,50000)
plot(prange, dbeta(prange, shape1=0.1+n[4]*0.5, shape2=0.1+n[4]*0.5), 
     type='l', col=mycolors[4], xlab="Taxon frequency (p)", 
     ylab="Probability(p|x,n)", main="Taxon frequency of 0.5")
lines(prange, dbeta(prange, shape1=0.1+n[3]*0.5, shape2=0.1+n[3]*0.5), col=mycolors[3])
lines(prange, dbeta(prange, shape1=0.1+n[2]*0.5, shape2=0.1+n[2]*0.5), col=mycolors[2])
lines(prange, dbeta(prange, shape1=0.1+n[1]*0.5, shape2=0.1+n[1]*0.5), col=mycolors[1])
legend(0.52, 170, legend=n, col=mycolors, lty=1, title="Read count", bty='n')

prange<-seq(0, 0.2, 0.001)
n<-c(100,1000,10000,50000)
plot(prange, dbeta(prange, shape1=0.1+n[4]*0.05, shape2=0.1+n[4]*0.95), 
     type='l', col=mycolors[4], xlab="Taxon frequency (p)", 
     ylab="Probability(p|x,n)", main="Taxon frequency of 0.05")
lines(prange, dbeta(prange, shape1=0.1+n[3]*0.05, shape2=0.1+n[3]*0.95), col=mycolors[3])
lines(prange, dbeta(prange, shape1=0.1+n[2]*0.05, shape2=0.1+n[2]*0.95), col=mycolors[2])
lines(prange, dbeta(prange, shape1=0.1+n[1]*0.05, shape2=0.1+n[1]*0.95), col=mycolors[1])
legend(0.1, 400, legend=n, col=mycolors, lty=1, title="Read count", bty='n')

n<-seq(100, 50000, 1000)
plot(n, qbeta(0.975, shape1=0.1+n*0.5, shape2=0.1+n*0.5) - 
       qbeta(0.025, shape1=0.1+n*0.5, shape2=0.1+n*0.5), xlab="Read count",
     ylab="Size of 95% credible interval for p", col=mycolors[1], type='l',
     ylim=c(0,0.2))
lines(n, qbeta(0.975, shape1=0.1+n*0.05, shape2=0.1+n*0.95) - 
       qbeta(0.025, shape1=0.1+n*0.05, shape2=0.1+n*0.95), col=mycolors[2])
legend(30000, 0.2, legend=c(0.5, 0.1), col=mycolors[1:2], lty=1, 
       title="Frequency", bty='n')
dev.print(pdf,"confidenceP.pdf", height=5,width=8)