code: def ref_fused_batch_normalization(
      x,
      beta,
      gamma,
      rmean,
      rvar,
      z,
      axes,
      decay_rate,
      eps,
      batch_stat,
      nonlinearity,
      output_stat,
  ):
      with nn.context_scope(cpu_context):
          xvar = nn.Variable.from_numpy_array(x)
          betavar = nn.Variable.from_numpy_array(beta)
          gammavar = nn.Variable.from_numpy_array(gamma)
          rmeanvar = nn.Variable.from_numpy_array(rmean)
          rvarvar = nn.Variable.from_numpy_array(rvar)
          if z is not None:
              zvar = nn.Variable.from_numpy_array(z)
          with nn.auto_forward():
              bn = F.batch_normalization(
                  xvar,
                  betavar,
                  gammavar,
                  rmeanvar,
                  rvarvar,
                  axes,
                  decay_rate,
                  eps,
                  batch_stat,
                  output_stat,
              )
              if z is None:
                  if output_stat:
                      y = bn[0]
                  else:
                      y = bn
              else:
                  if output_stat:
                      y = F.add2(bn[0], zvar)
                  else:
                      y = F.add2(bn, zvar)
              y = F.relu(y)
          rmean[:] = rmeanvar.d
          rvar[:] = rvarvar.d
          if output_stat:
              return y.d, bn[1].d, bn[2].d
          else:
              return y.d
  
file_path: nnabla/python/test/function/test_fused_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_fused_batch_normalization(x, beta, gamma, rmean, rvar, z, axes, decay_rate, eps, batch_stat, nonlinearity, output_stat) and returns: y,  bn,  bn, y