0

I have this custom function:

f <- function(N, zero_per, sed) {
  set.seed(sed)
  M <- N-(N*zero_per)
  y_g <- rgamma(M, shape=rnorm(M,1,2), rate=rnorm(M,1,2))
  df_g <- data.frame(y=y_g)
  df_l <- data.frame(y=rep(0, N-M))
  df <- rbind(df_g, df_l)
  return(df)
}

this function is part of a loop where I use a sequence for the second argument:

> sq <- seq(0.05, 0.7, 0.05)
> sq
 [1] 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70

the 12th position has value 0.6:

> sq[12]
[1] 0.6

However, when I call the function, once with the numeric argument 0.6, and once with sq[12], I get different results:

> head(f(5000, sq[12], 1))
          y
1       NaN
2       NaN
3       NaN
4 0.7547084
5 1.0643216
6       NaN
Warning message:
In rgamma(M, shape = rnorm(M, 1, 2), rate = rnorm(M, 1, 2)) : NAs produced

> head(f(5000, 0.6, 1))
          y
1       NaN
2       NaN
3       NaN
4 0.8990203
5 1.0719659
6       NaN
Warning message:
In rgamma(M, shape = rnorm(M, 1, 2), rate = rnorm(M, 1, 2)) : NAs produced

I use the same seed.

Strangely, this does not happen when I use 0.5:

> head(f(5000, sq[10], 1))
          y
1       NaN
2       NaN
3       NaN
4 0.4426849
5       NaN
6       NaN
Warning message:
In rgamma(M, shape = rnorm(M, 1, 2), rate = rnorm(M, 1, 2)) : NAs produced
> head(f(5000, 0.5, 1))
          y
1       NaN
2       NaN
3       NaN
4 0.4426849
5       NaN
6       NaN
Warning message:
In rgamma(M, shape = rnorm(M, 1, 2), rate = rnorm(M, 1, 2)) : NAs produced

The warnings can be ignored, they just happen because I stripped down my real function.

What is going on here?

spore234
  • 3,550
  • 6
  • 50
  • 76
  • 2
    Try `sq[12] - 0.6` and see that it is not zero. Classic example of [R FAQ 7.31](https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f), *"Why doesn’t R think these numbers are equal?"*. – r2evans Oct 17 '16 at 14:21

1 Answers1

1

The problem is that:

sq[12]==0.6
#[1] FALSE

While 

sq[10]==0.5
#[1] TRUE

To see that, do the following:

options(digits = 20)

sq[12]
#[1] 0.60000000000000009

sq[10]
#[1] 0.5
989
  • 12,579
  • 5
  • 31
  • 53