p <- c(construction = 0.10, retail = 0.05, agriculture = 0.20)
S <- expand.grid(c1 = 0:1, c2 = 0:1, c3 = 0:1) # 8 sample points
S$prob <- apply(S, 1, function(d) prod(ifelse(d == 1, p, 1 - p)))
S$Y <- rowSums(S[, 1:3]) # the function Y(E_i)
pY <- tapply(S$prob, S$Y, sum) # sum within {Y = y}
round(pY, 3); sum(pY) 0 1 2 3
0.684 0.283 0.032 0.001
[1] 1

