Round vector of numerics to integer while preserving their sum

Question

How to round floats to integers while preserving their sum? has the below answer written in pseudocode, which rounds a vector to integer values such that the sum of the elements in unchanged and the roundoff error is minimized. I'd like to implement this efficiently (i.e. vectorized if possible) in R.

For example, rounding these numbers yields a different total:

set.seed(1)
(v <- 10 * runif(4))
# [1] 2.655087 3.721239 5.728534 9.082078
(v <- c(v, 25 - sum(v)))
# [1] 2.655087 3.721239 5.728534 9.082078 3.813063
sum(v)
# [1] 25
sum(round(v))
# [1] 26

Copying pseudocode from answer for reference

// Temp array with same length as fn.
tempArr = Array(fn.length)

// Calculate the expected sum.
arraySum = sum(fn)

lowerSum = 0
-- Populate temp array.
for i = 1 to fn.lengthf
    tempArr[i] = { result: floor(fn[i]),              // Lower bound
                   difference: fn[i] - floor(fn[i]),  // Roundoff error
                   index: i }                         // Original index

    // Calculate the lower sum
    lowerSum = lowerSum + tempArr[i] + lowerBound
end for

// Sort the temp array on the roundoff error
sort(tempArr, "difference")

// Now arraySum - lowerSum gives us the difference between sums of these
// arrays. tempArr is ordered in such a way that the numbers closest to the
// next one are at the top.
difference = arraySum - lowerSum

// Add 1 to those most likely to round up to the next number so that
// the difference is nullified.
for i = (tempArr.length - difference + 1) to tempArr.length
    tempArr.result = tempArr.result + 1
end for

// Optionally sort the array based on the original index.
array(sort, "index")

Answer 1

Thanks for this useful function! Just to add to the answer, if rounding to the specified number of decimal places, the function can be modified:

smart.round <- function(x, digits = 0) {
  up <- 10 ^ digits
  x <- x * up
  y <- floor(x)
  indices <- tail(order(x-y), round(sum(x)) - sum(y))
  y[indices] <- y[indices] + 1
  y / up
}

Answer 2

In an even simpler form, I would say this algorithm is:

1. Start with everything rounded down
2. Round up the numbers with the highest fractional parts until the desired sum is reached.

This can be implemented in a vectorized way in R by:

1. Round down with floor
2. Order numbers by their fractional parts (using order)
3. Use tail to grab the indices of the elements with the k largest fractional parts, where k is the amount that we need to increase the sum to reach our target value
4. Increment the output value in each of these indices by 1

In code:

smart.round <- function(x) {
  y <- floor(x)
  indices <- tail(order(x-y), round(sum(x)) - sum(y))
  y[indices] <- y[indices] + 1
  y
}
v
# [1] 2.655087 3.721239 5.728534 9.082078 3.813063
sum(v)
# [1] 25
smart.round(v)
# [1] 2 4 6 9 4
sum(smart.round(v))
# [1] 25

Answer 3

Running total and diff based approach is much faster compared to smartRound by @josliber:

diffRound <- function(x) { 
  diff(c(0, round(cumsum(x)))) 
}

Here how results compare on 1m records (see details here: Running Rounding):

res <- microbenchmark(
  "diff(dww)" = x$diff.rounded <- diffRound(x$numbers) ,
  "smart(josliber)"= x$smart.rounded <- smartRound(x$numbers),
  times = 100
)

Unit: milliseconds
expr            min       lq        mean       median     uq       max       neval
diff(dww)       38.79636  59.70858  100.6581   95.4304    128.226  240.3088   100
smart(josliber) 466.06067 719.22723 966.6007   1106.2781  1177.523 1439.9360  100