Why is this Binning by Median code wrong?

Question

I was working on Binning by Mean, Median and Boundary in R.

# R CODE
a=c(20.5, 52.5, 62.6, 72.4, 104.8, 63.9, 35.3, 83.9, 37.4, 71.6, 74.6, 44.5, 66.6, 56.1, 45.3, 37.2)
a=sort(a)
binsize=4
median(a[1:4])
# BINS ARE
for(i in 1:length(a))
{
  if(i%%binsize==0)
  {
    print("HI")
    print(a[i-3])
    print(a[i-2])
    print(a[i-1])
    print(a[i])
  }
}

# BINNING BY MEAN
sum=0
for(i in 1:length(a))
{
  sum=sum+a[i]
    if(i%%binsize==0){
      avg=sum/binsize
      sum=0
      print(rep.int(avg,binsize))
    }
}

# BINNING BY MEDIAN
i=1
for(i in 1:length(a)){
  if(i%%binsize==0)
  {
    print(rep.int(median(a[i-3:i]),binsize))
  }
}

Can anyone tell me why binning by median is giving me wrong output. The median(a[i-3:i]) for 1st bin returns a value which is not same as median(a[1:4]) for 1st bin. Why?

Answer 1

The error is due to the well-known R gotcha that the : (colon operator, which calls seq()) takes higher precedence than arithmetic. Always parenthesize arguments to : if they involve arithmetic or are expressions: a[(i-3):i]

Your code a[i-3:i] doesn't do what you want it to do a[(i-3):i], it does a[i - (3:i)]). So the medians you are printing are for these slices:

4-3:1 # i.e. 1:3
8-3:1 # i.e. 5:7
12-3:1 # i.e. 9:11
16-4:1 # i.e. 13:15

PS some coding-style tips

• You don't need to iterate over all possible values of i and check them modulo binsize, just do:

  for(i in seq(binsize, length(a), by=binsize)) { ... }

so in your median case:

for(i in seq(binsize, length(a), by=binsize)) {
  print(rep.int(median(a[(i-3):i]),binsize))
}
[1] 36.25 36.25 36.25 36.25
[1] 48.9 48.9 48.9 48.9
[1] 65.25 65.25 65.25 65.25
[1] 79.25 79.25 79.25 79.25

• But in fact you can replace even that with:

  split(a, ceiling(seq_along(a)/binsize))

  as per the "Split a vector into chunks in R"

  To make it even clearer, you could define a helper function chunk <- function(x, binsize) { split(x, ceiling(seq_along(x)/binsize)) }

• Then you can replace the for-loop with sapply:

.

sapply(split(a, ceiling(seq_along(a)/binsize)), mean)
sapply(chunk(a,binsize), mean)
     1      2      3      4 
32.600 49.600 66.175 83.925

sapply(split(a, ceiling(seq_along(a)/binsize)), median)
sapply(chunk(a,binsize), median)
    1     2     3     4 
36.25 48.90 65.25 79.25

Much cleaner, easier to read, and prevents errors, right?

• The generalization of the colon operator is the seq() function, give it a read, it's pretty useful.