There is definitely reason for it. When you have a deep network of L layers with BN at almost every linear layer, and a large dataset of N samples. To get the accurate dataset statistics for each BN layer, you would need to activate the network O(NL) times, which is too much. Using running estimators to estimate the statistics of a dynamic (but converging) system is cheap and appears to give good results in the case of BN. So it is preferred in this scenario.