Bioequivalence and Bioavailability Forum

Hi mittyri,

❝ never trust in '==' or '!=' for floating point numbers comparison

❝ See for example here

Great essay!

  I have stopped reading Stephen King novels.
Now I just read C code instead.    Richard A. O’Keefe
  Being really good at C++ is like being really good
at using rocks to sharpen sticks.    Thant Tessman

In R-lingo:

test.while <- function(x0, x1, dx) {
  x    <- double()
  x[1] <- x0
  j    <- 1L
  while (max(x) < x1) {
    j <- j + 1L
    x[j] <- x[j-1] + dx
  }
  return(x)
}

test.repeat1 <- function(x0, x1, dx) {
  x <- x0
  repeat {
    if (max(x) >= x1) break
    x <- c(x, max(x) + dx)
  }
  return(x)
}

test.repeat2 <- function(x0, x1, dx) {
  x    <- double(round((x1-x0)/dx, 0) + 1)
  x[1] <- x0
  j    <- 1L
  repeat {
    j    <- j + 1L
    if (x[j-1] >= x1) break
    x[j] <- x[j-1] + dx
  }
  return(x)
}

test.for <- function(x0, x1, dx) {
  x <- double(ceiling((x1-x0)/dx) + 1)
  for (j in seq_along(x)) {
    x[j] <- x0 + (j-1)*dx
  }
  return(x)
}

test.seq <- function(x0, x1, dx) {
  x <- seq(from=x0, to=dx*ceiling((x1-x0)/dx), by=dx)
  return(x)
}

test.pretty <- function(x0, x1, dx) {
  n <- ceiling((x1-x0)/dx) + 1
  x <- pretty(x=c(x0, x1), n=n)
  return(x)
}

x0 <- 0
x1 <- 0.234567
dx <- 0.05

test.while(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25
test.repeat1(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25
test.repeat2(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25
test.for(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25
test.seq(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25
test.pretty(x0, x1, dx)
# [1] 0.00 0.05 0.10 0.15 0.20 0.25

Speed?

library(microbenchmark)
options(microbenchmark.unit="relative")
res <- microbenchmark(test.while(x0, x1, dx), test.repeat1(x0, x1, dx),
                      test.repeat2(x0, x1, dx), test.for(x0, x1, dx),
                      test.seq(x0, x1, dx), test.pretty(x0, x1, dx),
                      times=1000L, control=list("random", warmup=10))
print(res, digits=4)

Unit: relative
                     expr    min    lq  mean median  uq    max neval    cld
   test.while(x0, x1, dx)  2.333 2.249 2.233  2.222 2.1 0.9143  1000    d 
 test.repeat1(x0, x1, dx)  2.166 2.000 1.992  2.000 1.9 0.6714  1000   c   
 test.repeat2(x0, x1, dx)  1.333 1.375 1.342  1.333 1.3 0.3572  1000  b   
     test.for(x0, x1, dx)  1.000 1.000 1.000  1.000 1.0 1.0000  1000 a     
     test.seq(x0, x1, dx)  8.663 7.372 7.101  6.889 6.5 2.4857  1000     e
  test.pretty(x0, x1, dx) 10.662 8.746 8.393  8.111 7.6 3.5286  1000      f

Yes, we have a winner.
IMHO, seq() is most comfortable. One-liner: seq(0, 0.05*ceiling(0.234567/0.05), 0.05)
Good luck if you want nice ticks for a log-axis…

Hi Nobody and everybody else,

to illustrate it please have a look at this archive. Unzip to its own folder, then drag and drop the csv file onto the executable. You may need to give green light if your security package protests or something. There is no virus in it (at least there wasn't when the programmer uploaded it), and it does not phone home.

It demonstrates what the programmer is doing with a bootstrap solution for dissolution data. It provides 1.000.000 bootstraps of two datasets and illustrates the results in the html file.
The hot stuff, whose detail I will not discuss here (that's for now solely a matter between me and someone in a public institution who likes to look at dissolution and to reject dossiers in murderous ways ... I hope said person isn't reading this post), is the automatic generation of a series of graphs and plots and with an objective function, gamma. The user looks at these graphs, considers the gamma, and takes a decision. The user could in principle have a look at 1.000.000 such graphs and gammas (how about doing that in Prism, Nobody?).
The selection of curves to present is done in a pseudo-clever (actually pretty retarded, presently) fashion.

I hope this illustrates why a manual solution isn't viable (if not then then trust me, a manual solution is not viable. Not at all.)

For those who are into bootstraps, the result is at the end as a table, including raw bootstrap confidence intervals, Bias-corrected bootstrap interval, and bias-corrected and accelerated interval.

Hötzi: If you are reading this then I know you will absolutely torture my poor innocent little executable which is guilty of no fault . Yes, you can make it crash. Yes, the programmer (some wacko from Denmark) did not implement error trapping for whatever garbage the users fill in. No, this is not for the masses. Yes, the programmer has restricted the intended use to this dataset only.
By the way: The upload link function didn't work 2 minutes ago (it added bebac.at in the target string), so I pasted the raw link. Please forgive me.