code: def ref_batch_normalization_for_multiple_axes(
      x,
      beta,
      gamma,
      rmean,
      rvar,
      axes,
      decay_rate,
      eps,
      batch_stat,
      output_stat,
  ):
      axes = [a + x.ndim * (a < 0) for a in axes]
      reduc_axes = [i for i in range(x.ndim) if not i in axes]
      div_factor = 1
      for i in range(len(axes)):
          div_factor = div_factor * x.shape[axes[i]]
      reduc_axes = tuple(reduc_axes)
      m = x.size / div_factor
      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 6 for calling a global function with example usage: ref_batch_normalization_for_multiple_axes(x, beta, gamma, rmean, rvar, axes, decay_rate, eps, batch_stat, output_stat) and returns: y, y,  mean,  var