library(dplyr)
nSimulations <- 10000
expand.grid(sim = 1:nSimulations) %>%
group_by(sim) %>%
mutate(x = runif(1), ok1 = x < 0.5,
y = runif(1,min=0, max =1-x), ok2 = y < 0.5,
z= 1 - x -y, ok3 = z < 0.5,
ok = ok1 && ok2 && ok3) %>%
ungroup() %>%
summarise(prob = sum(ok) / nSimulations)
## Source: local data frame [1 x 1]
##
## prob
## 1 0.1941