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…