## This is free and unencumbered software released into the public domain.
##
## Anyone is free to copy, modify, publish, use, compile, sell, or
## distribute this software, either in source code form or as a compiled
## binary, for any purpose, commercial or non-commercial, and by any
## means.
##
## In jurisdictions that recognize copyright laws, the author or authors
## of this software dedicate any and all copyright interest in the
## software to the public domain. We make this dedication for the benefit
## of the public at large and to the detriment of our heirs and
## successors. We intend this dedication to be an overt act of
## relinquishment in perpetuity of all present and future rights to this
## software under copyright law.
##
## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
## EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
## MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
## IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
## OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
## ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
## OTHER DEALINGS IN THE SOFTWARE.
##
## For more information, please refer to <http://unlicense.org/>

badGenerator <- function(mx) {
  r <- 2^ceiling(log(mx) / log(2))
  sample(1:r, 1) %%  mx + 1
}

badSample <- function(n) {
    acc <- c()
    for( ii in 1:length(n)){
      index <- badGenerator(length(n))
      acc <- c(acc, n[index])
      n <- n[-index]
    }
    acc
}

maxVal <- 128
bias_pos_count <- rep(0, times=maxVal)
unbias_pos_count <- rep(0, times=maxVal)
trials <- 100000

for( ii in 1:trials ){
    bias_pos_count <- bias_pos_count +  badSample( 1:maxVal )
    unbias_pos_count <- unbias_pos_count + sample( 1:maxVal )
}

bias_means <- bias_pos_count / trials
unbias_means <- unbias_pos_count / trials

shift <- bias_means - unbias_means
dev <- sd(shift)
deviation <- shift / dev


color_map = heat.colors(maxVal)[rank(shift)]

#barplot(bias_means, xlab="Index", ylab="Deviation", bty="n", main="Mean shift (trials = 10^5)", col = color_map)
#barplot(deviation, xlab="Index", ylab="Deviation", bty="n", main="Mean shift (trials = 10^5)", col = color_map)
barplot(deviation, xlab="Index", ylab="Deviation", bty="n", main="Mean shift from equal likelihood (trials = 10^5)", col = color_map)