Search code examples
rroundingmarkov-chains

Rounding Stochastic Variables in R


When converting decimal values in a transition matrix to fractions with a specific denominator, how can I insure the resulting rows sum to 1? I understand that this is an issue relating to the choice of denominator, in this case x/36; however, I want to force rounding up of e.g. 0/36 to 1/36 to insure the resulting rows sum to 1, even if the error is large. Ideally the forced rounding would correspond to the highest likelihood from the transition matrix. The application is related to manual game design in which the transition matrix is used to represent states in the environment and is traversed by the sum of two six-sided dice with 36 possible combinations.

I have tried to use different rounding algorithms on the decimal values including the base round() and smart.round as documented here Round vector of numerics to integer while preserving their sum, but I am largely at a loss for how to modify smart.round to fit my particular use case. Am also using function documented here Function in R to convert a decimal to a fraction with a specified denominator for the conversion of decimals to fractions.

# transition matrix

transmat <- structure(c(0.730926989335521, 0.474022677271958, 0.0091743119266055, 
                        0.326860841423948, 0.97411003236246, 0.926605504587156, 0.926605504587156, 
                        0.961722488038278, 0.926605504587156, 0.238720262510254, 0.40632932814351, 
                        0.926605504587156, 0.326860841423948, 0.00323624595469256, 0.0091743119266055, 
                        0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                        7.45767767917071e-05, 0.0170925706549332, 0.0091743119266055, 
                        0.00323624595469256, 0.00323624595469256, 0.0091743119266055, 
                        0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                        0.00753225445596242, 0.0340159079370452, 0.0091743119266055, 
                        0.00323624595469256, 0.00323624595469256, 0.0091743119266055, 
                        0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                        0.00753225445596242, 0.0340159079370452, 0.0091743119266055, 
                        0.00323624595469256, 0.00323624595469256, 0.0091743119266055, 
                        0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                        7.45767767917071e-05, 0.0170925706549332, 0.0091743119266055, 
                        0.00323624595469256, 0.00323624595469256, 0.0091743119266055, 
                        0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                        7.45767767917071e-05, 0.0170925706549332, 0.0091743119266055, 
                        0.00323624595469256, 0.00323624595469256, 0.0091743119266055, 
                        0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                        0.00753225445596242, 0.00016923337282112, 0.0091743119266055, 
                        0.326860841423948, 0.00323624595469256, 0.0091743119266055, 0.0091743119266055, 
                        0.00478468899521532, 0.0091743119266055, 0.00753225445596242, 
                        0.00016923337282112, 0.0091743119266055, 0.00323624595469256, 
                        0.00323624595469256, 0.0091743119266055, 0.0091743119266055, 
                        0.00478468899521532, 0.0091743119266055), dim = c(9L, 9L), dimnames = list(
                          c("00", "02", "13", "18", "25", "50", "60", "80", "95"), 
                          c("00", "02", "13", "18", "25", "50", "60", "80", "95")))

# smart.round and function for decimal to fractions conversion with form x/36

smart.round <- function(x) {
  y <- floor(x)
  indices <- tail(order(x-y), round(sum(x)) - sum(y))
  y[indices] <- y[indices] + 1
  y
}

f = function(x) {paste0(smart.round(x * 36), "/", 36)}

# convert transmat values to fractions

fractions <- apply(transmat, c(1,2), f)

# > sessionInfo()
# R version 4.2.2 (2022-10-31 ucrt)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 10 x64 (build 19045)
# 
# Matrix products: default
# 
# locale:
#   [1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
# [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
# [5] LC_TIME=English_United States.utf8    
# 
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# loaded via a namespace (and not attached):
#   [1] fansi_1.0.3      utf8_1.2.2       R6_2.5.1         lifecycle_1.0.3  magrittr_2.0.3  
# [6] scales_1.2.1     pillar_1.8.1     rlang_1.0.6      cli_3.4.1        rstudioapi_0.14 
# [11] paletteer_1.5.0  vctrs_0.5.0      rematch2_2.1.2   tools_4.2.2      glue_1.6.2      
# [16] munsell_0.5.0    compiler_4.2.2   pkgconfig_2.0.3  colorspace_2.0-3 tibble_3.1.8


Solution

  • The R package pomdp implements smart.round in the round_stochastic function, which rounds the values of a matrix row such that the row sum is preserved.

    Note also that in my original question I was not applying the function f correctly to the transition matrix. You must transpose the matrix with the call to apply and use the margin option correctly. Thus,

    fractions <- apply(transmat, c(1,2), f)
    

    becomes:

    fractions <- t(apply(transmat, c(1), f))
    
        library(pomdp)
        
        # transition matrix
        
        transmat <- structure(c(0.730926989335521, 0.474022677271958, 0.0091743119266055, 
                                0.326860841423948, 0.97411003236246, 0.926605504587156, 0.926605504587156, 
                                0.961722488038278, 0.926605504587156, 0.238720262510254, 0.40632932814351, 
                                0.926605504587156, 0.326860841423948, 0.00323624595469256, 0.0091743119266055, 
                                0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                                7.45767767917071e-05, 0.0170925706549332, 0.0091743119266055, 
                                0.00323624595469256, 0.00323624595469256, 0.0091743119266055, 
                                0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                                0.00753225445596242, 0.0340159079370452, 0.0091743119266055, 
                                0.00323624595469256, 0.00323624595469256, 0.0091743119266055, 
                                0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                                0.00753225445596242, 0.0340159079370452, 0.0091743119266055, 
                                0.00323624595469256, 0.00323624595469256, 0.0091743119266055, 
                                0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                                7.45767767917071e-05, 0.0170925706549332, 0.0091743119266055, 
                                0.00323624595469256, 0.00323624595469256, 0.0091743119266055, 
                                0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                                7.45767767917071e-05, 0.0170925706549332, 0.0091743119266055, 
                                0.00323624595469256, 0.00323624595469256, 0.0091743119266055, 
                                0.0091743119266055, 0.00478468899521532, 0.0091743119266055, 
                                0.00753225445596242, 0.00016923337282112, 0.0091743119266055, 
                                0.326860841423948, 0.00323624595469256, 0.0091743119266055, 0.0091743119266055, 
                                0.00478468899521532, 0.0091743119266055, 0.00753225445596242, 
                                0.00016923337282112, 0.0091743119266055, 0.00323624595469256, 
                                0.00323624595469256, 0.0091743119266055, 0.0091743119266055, 
                                0.00478468899521532, 0.0091743119266055), dim = c(9L, 9L), dimnames = list(
                                  c("00", "02", "13", "18", "25", "50", "60", "80", "95"), 
                                  c("00", "02", "13", "18", "25", "50", "60", "80", "95")))
        
        # use round_stochastic function (pomdp) to round rows. Fraction conversion to x/36 form.
        
        f = function(x) {paste0(round_stochastic(x * 36, digits = 0), "/", 36)}
        
        # convert transmat values to fractions
        
        fractions <- t(apply(transmat, c(1), f))