rmbeta<-
function (n = 1000, mean = 0.5, sd = 0.2)
{
a <- mean/(sd * sd)
b <- (1 - mean)/(sd * sd)
rbeta(n, a, b)
}
masks_work<-rmbeta(10000,0.5,0.5)
hist(masks_work)
masks_worn_correctly<-rmbeta(1000,0.5,0.5)
masking<-masks_work*masks_worn_correctly
hist(masking)
ventillation<-rmbeta(10000,0.9,0.2)
hist(ventillation)
distancing<-rmbeta(10000,0.5,1)
hist(distancing)
d<-data.frame(masking,ventillation,distancing)
hist(apply(d,1,max))