pop1 = matrix(rbinom(10^8,size=1,p=0.5), ncol=100000)
pop2 = matrix(rbinom(10^8,size=1,p=0.5), ncol=100000)

p1 = apply (pop1, 2, mean)
q1 = 1-p1
p2 = apply (pop2, 2, mean)
q2=1-p2


sc1.1 = pop1 %*% log(p1) + (1-pop1) %*% log(q1)
sc1.2 = pop1 %*% log(p2) + (1-pop1) %*% log(q2)
sc2.1 = pop2 %*% log(p1) + (1-pop2) %*% log(q1)
sc2.2 = pop2 %*% log(p2) + (1-pop2) %*% log(q2)

# make sure you set the local directory to find the pdf
pdf("privacy.pdf")
par(mfrow = c(2,2))
plot (density(sc1.1),xlim = c(min(c(sc1.2,sc2.1)), max(c(sc1.1,sc2.2))), main = "case log likelihood")
lines (density(sc1.2),col="red")
legend("topright",legend=c("cases", "controls"), lty = c(1,1), col=c(1,"red"))

plot (density(sc2.2),xlim = c(min(sc2.1), max(sc2.2)), main = "control log likelihood")
lines (density(sc2.1),col="red")
legend("topright",legend=c("controls", "cases" ), lty = c(1,1), col=c(1,"red"))

none = matrix(rbinom(10^7,size=1,p=0.5), ncol=100000)

none.1 = none %*% log(p1) + (1-none) %*% log(q1)
none.2 = none %*% log(p2) + (1-none) %*% log(q2)

plot (density(none.1),xlim = c(min(c(sc1.2,sc2.1)), max(c(sc1.1,sc2.2))),col="blue", main = "new sample log likelihood")
lines (density(none.2),col="green")
legend("topright",legend=c("controls", "cases" ), lty = c(1,1), col=c("blue","green"))
dev.off()