alpha <- 2.92 # From O'Quigley(2024)#
beta <- 16.23 # From O'Quigley(2024)
library(MASS)
set.seed(123456)
p <- (0:1000)/1000
rbetas <- rbeta(38, alpha, beta)
truehist(rbetas, h = 0.05, xlim = c(0, 1), ylim = c(0, 7), prob = TRUE,
main = "Histogram of a sample of size 38 \n from O'Quigley's fitted beta distribution",
sub = "alpha = 2.92, beta = 16.23",
xlab = "Nurse specific probability of shift rostering", ylab = "Probability density")
lines(p, dbeta(p, alpha, beta), col = "red", lwd = 2)
lines(x = c(0, 1), y = c(0, 0)); lines(x = c(0, 0), y = c(0, 7))
