Introduction

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))