A Monte Carlo simulator is developed to investigate the low resistance state (LRS) retention of HfOx based resistive switching memory. The simulator tracks the evolution of oxygen ions and oxygen vacancies by choosing the occurrence of the generation/recombination/migration processes based on their corresponding probability at each simulation step. The current through the resulting conductive filament (CF) configuration is then calculated by a trap-assisted-tunneling solver. The simulation results are corroborated with experimental data. The LRS retention is found to be CF size dependent, and its variability is suggested to be intrinsic due to the stochastic nature of the CF dissolution. The tail bits of the failure time distribution become a limiting factor of the device reliability in a large memory array.