I had an odd problem where I wanted to round numbers probabilistically. As an example if we have the number 7. 5 we should round to 7 half the time and to 8 half the time. This tie breaking method is commonly known as stochastic rounding. In contrast to this practice I want to round a number like 7.25 probabilistically as well. For instance, I would round to 7.25 to 7 three quarters of the time and to 8 a quarter of the time. I couldn’t find anything quite like this in R so I wrote the small function below that will handle positive and negative numbers.
StochRound = function(x){
## extract the decimal portion
q = abs(x - trunc(x))
## draw a value 0 or 1 with probability
## based on how close we already are
adj = sample(0:1, size = 1, prob = c(1 - q, q))
## make it negative if x is
if(x < 0) adj = adj * -1
## return our new value
trunc(x) + adj
}
This seems to do the trick. To check and make sure that it doesn't create a bias of some sort we can repeat this process several 1,000 times and make sure that our mean isn't changing.
## extract the decimal portion
q = abs(x - trunc(x))
## draw a value 0 or 1 with probability
## based on how close we already are
adj = sample(0:1, size = 1, prob = c(1 - q, q))
## make it negative if x is
if(x < 0) adj = adj * -1
## return our new value
trunc(x) + adj
}
This seems to do the trick. To check and make sure that it doesn't create a bias of some sort we can repeat this process several 1,000 times and make sure that our mean isn't changing.
x <- mean="1)</font" n="100000," rnorm="">
mean(x)
x <- font="" sapply="" stochround="" x="">
mean(x)
mean(x)
x <- font="" sapply="" stochround="" x="">
mean(x)
mean before was 1.002174 and after rounding we were at 1.00278 well within what we expect. Finally we might even just look at the distribution before and after rounding:
I think that is good enough to go ahead and use it.
The R code is really nice and simple, cool. And I really like the verification approach.
ReplyDelete