shift_length <- 480 #minutes

N <- 1e4 #Number of simulations

#Concentrations and durations for each scenario

c1 <- rlnorm(N, log(2.3), log(1.5))

t1 <- rnorm(N, 30, 5)

c2 <- rlnorm(N, log(0.8), log(2.2))

t2 <- rnorm(N, 120, 20)

c3 <- rlnorm(N, log(1.5), log(1.7))

t3 <- rnorm(N, 60, 10)

#Proportion of shift covered by task samples

(mean(t1) + mean(t2) + mean(t3)) / shift_length

#TWA Estimation

TWA <- (c1 * t1 + c2 * t2 + c3 * t3) / shift_length

#Visualisation

plot(TWA,

main = "Simulated TWAs",

ylab = "Concentration (mg/m3)")

abline(1,0, col="red", lwd = 2)

hist(TWA,

breaks = 30,

main = "Estimated TWAs",

xlab = "Concentration (mg/m3)",

ylab = "Likelihood",

col = "lightblue",

border = "black",

probability = TRUE)

abline(v=1, col="red", lwd = 2)

#Summarising simultaed TWA

print(mean(TWA)) #Expected TWA

print(quantile(TWA, probs = seq(0, 1, 0.25))) #TWA ranges

print(mean(TWA >=1)) # Probability of exceedance

#Concentration Hours of Each Task

#Bigger numbers means heigher priority for control

mean(c1*t1)

mean(c2*t2)

mean(c3*t3)