code: def double_backward_for_batch(
      g_dx0, g_db0, g_dg0, dy, x0, b0, g0, rm, rv, axes, decay_rate, eps
  ):
      # Prerequisite
      # axes reduced and denominator
      axes0 = [a for a in range(x0.ndim)]
      axes = [a + x0.ndim * (a < 0) for a in axes]
      axes = list(set(axes0) - set(axes))
      F_sum = partial(F.sum, axis=axes, keepdims=True)
      F_mean = partial(F.mean, axis=axes, keepdims=True)
      de = np.prod([x0.shape[a] for a in axes])  # denominator
      bm = F_mean(x0)  # batch mean
      bv = F_mean(x0**2.0) - bm**2.0  # batch variance
      x0_bm = x0 - bm  # x0 - batch mean
      v_eps_r1 = (bv + eps) ** -1.0  # (variance + eps) * (-1)
      v_eps_rsqrt1 = (bv + eps) ** (
          -1.0 / 2.0
      )  # (variance + eps) * (-1/2)
      v_eps_rsqrt3 = v_eps_rsqrt1**3.0  # (variance + eps) * (-3/2)
  
      # common factors
      dy_x0_bm_sum = F_sum(dy * x0_bm)
      dy_sum = F_sum(dy)
      g_dx0_x0_bm_sum = F_sum(g_dx0 * x0_bm)
      g_dx0_sum = F_sum(g_dx0)
  
      # wrt. x
      # from dx
      dv = (-1.0 / 2.0) * g0 * v_eps_rsqrt3 * F_sum(dy * x0_bm)
      g_dx0_dy_sum = F_sum(g_dx0 * dy)
  
      g1 = (-1.0 / de) * v_eps_rsqrt3 * g_dx0_dy_sum * g0 * x0_bm
      g2 = (
          (1.0 / de)
          * g0
          * g_dx0_x0_bm_sum
          * v_eps_rsqrt3
          * (
              1.0 / de * dy_sum
              - dy
              + (3.0 / de) * v_eps_r1 * dy_x0_bm_sum * x0_bm
          )
      )
      g2 += (2.0 / de) * dv * (g_dx0 - (1.0 / de) * g_dx0_sum)
      g3 = (
          (1.0 / de**2.0)
          * g_dx0_sum
          * dy_sum
          * g0
          * v_eps_rsqrt3
          * x0_bm
      )
      g_x0 = g1 + g2 + g3
  
      # from gamma
      t1 = (dy - dy_sum / de) * v_eps_rsqrt1
      t2 = (-1.0 / de) * dy_x0_bm_sum * v_eps_rsqrt3 * x0_bm
      g_x0 += g_dg0 * (t1 + t2)
  
      # wrt. beta
      # zero, do nothing
  
      # wrt. gamma
      t1 = dy * v_eps_rsqrt1
      t2 = (-1.0 / de) * dy_x0_bm_sum * v_eps_rsqrt3 * x0_bm
      t3 = (-1.0 / de) * dy_sum * v_eps_rsqrt1
      g_g0 = F_sum(g_dx0 * (t1 + t2 + t3))
  
      # wrt. dy
      t1 = g_dx0 * g0 * v_eps_rsqrt1
      t2 = -(1.0 / de) * g0 * v_eps_rsqrt3 * g_dx0_x0_bm_sum * x0_bm
      t3 = -(1.0 / de) * g0 * v_eps_rsqrt1 * g_dx0_sum
      x0_hat = x0_bm * v_eps_rsqrt1
      g_dy = (t1 + t2 + t3) + g_dg0 * x0_hat + g_db0
  
      return g_dy, g_x0, None, g_g0
  
imports: from functools import partial
  import nnabla
  import nnabla.function
  import nnabla.functions
  import numpy
  from nnabla.function import PythonFunction
  from .utils import force_tuple
label_desc: Use image node_id 1 for calling a global function with example usage: double_backward_for_batch(g_dx0, g_db0, g_dg0, dy, x0, b0, g0, rm, rv, axes, decay_rate, eps) and returns: g_dy,  g_x0,  None,  g_g0