code: def ref_batch_normalization(
      x,
      beta,
      gamma,
      rmean,
      rvar,
      axes,
      decay_rate,
      eps,
      batch_stat,
      output_stat,
  ):
      beta, gamma, rmean, rvar
  
      assert len(axes) == 1
      reduc_axes = list(range(x.ndim))
      del reduc_axes[axes[0]]
      reduc_axes = tuple(reduc_axes)
      m = x.size / x.shape[axes[0]]
      if batch_stat:
          mean = x.mean(reduc_axes, keepdims=True)
          var = x.var(reduc_axes, keepdims=True)
          rmean[...] = decay_rate * rmean + (1 - decay_rate) * mean
          rvar[...] = decay_rate * rvar + (1 - decay_rate) * var * m / (
              m - 1
          )
      else:
          mean = rmean
          var = rvar
      x_hat = (x - mean) / np.sqrt(var + eps)
      y = x_hat * gamma + beta
      if output_stat:
          return y, mean, var
      return y
  
file_path: nnabla/python/test/function/test_batch_normalization.py
imports: from six.moves import range
  import pytest
  import numpy
  import nnabla
  import nnabla.functions
  from nbla_test_utils import list_context
  from nnabla.testing import assert_allclose
label_desc: Use image node_id 1 for calling a global function with example usage: ref_batch_normalization(x, beta, gamma, rmean, rvar, axes, decay_rate, eps, batch_stat, output_stat) and returns: y, y,  mean,  var