How to simulate the distribution of a diagnostic test I have a parameter described as mean (95% CI low-hi) and would like to build a distribution that will have approximately that mean and lower/upper bounds. The constraints are that these are usually non-normal and all values have to be between 0 and 1. How would I go about that? 
One brute-force method I thought about is below, but I was wondering if there is a proper or, at least, more elegant way to do this. 
Note: A previous version of this question mentioned that these parameters are diagnostic test characteristics such as sensitivity and specificity. 
bMoms <- function(m, sd) {
  ## the moments of the beta distribution   
  v <- sd^2
  t <- (1 - m)/v
  a <- m^2 * t
  b <- m * t - a
  ## method of moments
  t. <- m * (1 - m)/v - 1
  a. <- m * t.
  b. <- (1 - m) * t.
  return(c(a, b, a., b.))
}
m = (lo + 4 * m + hi)/6
m. <- (m - lo)/(hi - lo)  #normalize m to scale 0-1
sd. <- 1/6
ab <- bMoms(m., sd.)  # the moments of the beta distribution
df. <- data.frame(x = rbeta(n, ab[1], ab[2]))
df. <- within(df., x <- lo + x * (hi - lo))

 A: Simulation is a viable means of inspecting the operating characteristics of a diagnostic test. While preliminary studies present sensitivity and specificity with 95% confidence intervals, these CIs do not actually represent a Bayesian probability in what such values are. Thus when you set out to conduct your simulation study, it is better to consider the test characteristics as fixed or given, and then to simulate a random process of the disease status and test evaluation according to those parameters. By doing this, you can appropriately summarize the test performance (and its uncertainty) in a variety of scenarios, such as an elevated prevalence of disease or in larger samples. The basic parameters are given as follows:
\begin{array}{c|ccc}
 & D & \bar{D} & \\ \hline
T & TP & FP& n_t\\
\bar{T} & FN & TN & n-n_t\\
 & n_d & n-n_d
\end{array}
As an example, suppose a test was developed and studied in a cohort with 40% sensitivity ($TP/n_d$) and 90% specificity ($TN/(n-n_d)$). Suppose further I consult the literature and find that the originating sample was 100 patients with a disease prevalence of 10% (say, in a hospital setting). Then, the table above boils down to 
\begin{array}{c|ccc}
 & D & \bar{D} & \\ \hline
T & 4 & 9 & 13\\
\bar{T} & 6 & 81 & 87\\
 & 10 & 90
\end{array}
But if we wish to simulate these outcomes, the $TN$ and the $FN$ are merely binomial probability models. Using R I can simulate this contingency table with:
sim <- function(n, pd, se, sp) {
  nd <- rbinom(1, n, pd)
  tp <- rbinom(1, nd, se)
  tn <- rbinom(1, n-nd, sp)
  matrix(c(tp, nd-tp, n-nd-tn, tn), 2, 2)
}
set.seed(123)
sim(n=100, pd=0.1, se=0.4, sp=0.9)

gets:
> sim(n=100, pd=0.1, se=0.4, sp=0.9)
     [,1] [,2]
[1,]    4    8
[2,]    4   84

which shows some stochastic variability from the design. Performing many 1000s of replications:
OUT <- replicate(1000, sim(n=100, pd=0.1, se=0.4, sp=0.9))
sens <- OUT[1,1,] / colSums(OUT[,1,])
mean(sens) + c('lower 95%'=-1, 'mean'=0, 'upper 95%'=1) %o%  qnorm(0.975) * sd(sens)

Which gives 
> mean(sens) + c('lower 95%'=-1, 'mean'=0, 'upper 95%'=1) %o%  qnorm(0.975) * sd(sens) 
                [,1]
lower 95% 0.06350984
mean      0.40045478
upper 95% 0.73739973

So as expected this design is inefficient for evaluating test-reliability. We did not generate variability in the test, but rather simulated the test outcome as a random process which is accurately summarized in the sensitivity and specificity and varies under experimentation a function of experimental design.
