3

A bit of a bizarre case: I have code that creates 2x2 matrices and populates them with particular values, then generates other matrices with lengths equal to the cell values in the first matrices. When I look at the values in the matrix, I get one result - and when I go to generate the lists, they act as a completely different number.

For example, the relevant code, running on R 3.2.2

N <-100 
pa <-.3
pna <- .4
pb <-.1
pnb <-.2

table.nH <- matrix(nrow=2,ncol=2)

table.nH[1,1] <- (1-pna)*(1-pnb)*N

table.nH[2,1] <- (1-pa)*(1-pnb)*N

table.nH[1,1] #This should be 48

> [1] 48

table.nH[2,1] #This should be 56

> [1] 56

pass.list <- matrix(rep(c(0,0,0,0),table.nH[1,1]),ncol=4,byrow=TRUE)

fail.list <- matrix(rep(c(0,0,0,0),table.nH[2,1]),ncol=4,byrow=TRUE)

length(pass.list[,1]) #Number of rows *should be* 48

> [1] 48

length(fail.list[,1]) #Number of rows *should be* 56

> [1] 55

I can't for the life of me figure out why one works but the other doesn't. It works if I hard-code the 56 (i.e. rep(c(0,0,0,0),56)) but not if I reference the 56 in the matrix. The actual table in question has no NA values, and this isn't unique (I see it on other versions of the same problem for some reason). However, it's also not universal: some other tables generated this way work exactly as expected.

What's going on, and how can I fix it?

Josh O'Brien
  • 159,210
  • 26
  • 366
  • 455
Brian D
  • 33
  • 3

1 Answers1

3

The reason for this lies in the behaviour of rep. In the documentation, it says

Non-integer values of times will be truncated towards zero.

Now, if you look at the values in your matrix in detail, you will note, that one of the values is not an integer:

format(table.nH, digits = 18)
#      [,1]                  [,2]                 
# [1,] "48.0000000000000000" "                 NA"
# [2,] "55.9999999999999929" "                 NA"

If table.nH[2, 1] is rounded towards zero, you will get 55.

You can solve this by rounding to an integer

table.nH <- round(table.nH, 0)
format(table.nH, digits = 18)
#      [,1] [,2]
# [1,] "48" "NA"
# [2,] "56" "NA"

In case you wonder why table.nH[2, 1] is not an integer: this has to do with the way floating point numbers are represented. There is a question on this topic that has an excellent answer.

Community
  • 1
  • 1
Stibu
  • 15,166
  • 6
  • 57
  • 71
  • Thank you very kindly; your answer explained not only what was wrong, but why I overlooked it, and the answer you linked makes perfect sense (...in hindsight). – Brian D Dec 04 '15 at 04:20
  • You are welcome. If the answer helped, please cosider accepting it. – Stibu Dec 04 '15 at 06:02