# Helper function to calculate the F statistic for two independent samples.
# (In fact the code calculates the t statistic and then squares it.)
F.stat <- function (m1, sd1, n1, m2, sd2, n2) {
  var.p <- (((sd1 ^ 2) * (n1 - 1)) + ((sd2 ^ 2) * (n2 - 1))) / (n1 + n2 - 2)
  se.p <- sqrt((var.p / n1) + (var.p / n2))
  result <- ((m1 - m2) / se.p) ^ 2

  return(result)
}

# dp: number of decimal places
# nsamp: number of samples to try

perfectMSD <- function (m1, sd1, n1, m2, sd2, n2, dp=2, nsamp=1000000) {
  p10 <- 10 ^ dp
  m1.r <- round(m1, dp)
  sd1.r <- round(sd1, dp)
  m2.r <- round(m2, dp)
  sd2.r <- round(sd2, dp)

# F.pub is the F statistic that we will get from the rounded means and sds.
  F.pub <- round(F.stat(m1.r, sd1.r, n1, m2.r, sd2.r, n2), dp)

# Calculate the minimum and maximum F statistics.
# For simplicity of the code with min/max and sign changes, we use 0 as one of the means.
  md <- (m1.r - m2.r)
  d <- 0.49999999 / (10 ^ dp)
  dm <- 2 * d * sign(md)
  md.min <- md - dm
  F.min <- round(F.stat(md.min, sd1.r + d, n1, 0, sd2.r + d, n2), dp)
  md.max <- md + dm
  F.max <- round(F.stat(md.max, sd1.r - d, n1, 0, sd2.r - d, n2), dp)

# Generate plausible unrounded mean and SD values.
  m1.ns <- rep(m1.r, nsamp) + ((runif(nsamp) - 0.5) / p10)
  sd1.ns <- rep(sd1.r, nsamp) + ((runif(nsamp) - 0.5) / p10)
  m2.ns <- rep(m2.r, nsamp) + ((runif(nsamp) - 0.5) / p10)
  sd2.ns <- rep(sd2.r, nsamp) + ((runif(nsamp) - 0.5) / p10)

# Sanity check that generated numbers round to published ones.
  m1.sc <- sum(round(m1.ns, dp) - m1.r)
  sd1.sc <- sum(round(sd1.ns, dp) - sd1.r)
  m2.sc <- sum(round(m2.ns, dp) - m2.r)
  sd2.sc <- sum(round(sd2.ns, dp) - sd2.r)
  if ((m1.sc != 0) || (sd1.sc != 0) || (m2.sc != 0) || (sd2.sc != 0)) {
    stop("Bad rounding of mean or SD")
  }

# Generate the F values from our "ragged" means and SDs.
  F.ns <- F.stat(m1.ns, sd1.ns, n1, m2.ns, sd2.ns, n2)

# Round the F values and see how many match the published one.
  F.round <- round(F.ns, dp)
  F.minobs <- min(F.round)
  F.maxobs <- max(F.round)
  F.exact <- (F.round == F.pub)
  frac.exact <- sum(F.exact) / nsamp
  return(list(F.pub, F.min, F.max, frac.exact, F.minobs, F.maxobs))
}

dp <- 3
nsamp <- 1000000
set.seed(1)

catFrac <- function (m1, sd1, n1, m2, sd2, n2, dp=2, nsamp=1000000) {
  perf <- perfectMSD(m1, sd1, n1, m2, sd2, n2, dp, nsamp)
  format <- paste("%.", dp, "f", sep="")
  cat("DP=", dp, " n1=", n1, " n2=", n2,
      " FPerf=", sprintf(format, (perf[[1]])),
      " FMin=", sprintf(format, (perf[[2]])),
      " FMax=", sprintf(format, (perf[[3]])),
      " Matched=", sprintf(format, perf[[4]] * 100), "%",
      " FMinObs=", sprintf(format, (perf[[5]])),
      " FMaxObs=", sprintf(format, (perf[[6]])),
      "\n", sep="")
}

catFrac(8.10, 1.20, 50, 7.28, 1.46, 50, 2, nsamp)
catFrac(6.84, 1.46, 50, 5.20, 1.92, 50, 2, nsamp)
catFrac(7.30, 1.46, 50, 6.44, 1.68, 50, 2, nsamp)
catFrac(6.56, 1.72, 50, 5.88, 2.19, 50, 2, nsamp)
catFrac(7.46, 1.67, 50, 6.48, 2.11, 50, 2, nsamp)