Solved by mufeili in post #2 
              
              
                I think I have figured it out what’s going on. According to backward function for multinomial, they used -r/(1e-6+p) directly in the backward process.