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