code: def reml(
      sigma,
      components,
      design=None,
      n=1,
      niter=128,
      penalty_cov=np.exp(-32),
      penalty_mean=0,
  ):
      """
  
      Adapted from spm_reml.m
  
      ReML estimation of covariance components from sigma using design matrix.
  
      INPUTS:
          sigma        -- m-by-m covariance matrix
          components   -- q-by-m-by-m array of variance components
                          mean of sigma is modeled as a some over components[i]
          design       -- m-by-p design matrix whose effect is to be removed for
                          ReML. If None, no effect removed (???)
          n            -- degrees of freedom of sigma
          penalty_cov  -- quadratic penalty to be applied in Fisher algorithm.
                          If the value is a float, f, the penalty is
                          f * identity(m). If the value is a 1d array, this is
                          the diagonal of the penalty.
          penalty_mean -- mean of quadratic penalty to be applied in Fisher
                          algorithm. If the value is a float, f, the location
                          is f * np.ones(m).
  
  
      OUTPUTS:
          C            -- estimated mean of sigma
          h            -- array of length q representing coefficients
                          of variance components
          cov_h        -- estimated covariance matrix of h
      """
  
      # initialise coefficient, gradient, Hessian
  
      Q = components
      PQ = np.zeros(Q.shape)
  
      q = Q.shape[0]
      m = Q.shape[1]
  
      # coefficient
      h = np.array([np.diag(Q[i]).mean() for i in range(q)])
  
      ## SPM initialization
      ## h = np.array([np.any(np.diag(Q[i])) for i in range(q)]).astype(np.float64)
  
      C = np.sum([h[i] * Q[i] for i in range(Q.shape[0])], axis=0)
  
      # gradient in Fisher algorithm
  
      dFdh = np.zeros(q)
  
      # Hessian in Fisher algorithm
      dFdhh = np.zeros((q, q))
  
      # penalty terms
  
      penalty_cov = np.asarray(penalty_cov)
      if penalty_cov.shape == ():
          penalty_cov = penalty_cov * np.identity(q)
      elif penalty_cov.shape == (q,):
          penalty_cov = np.diag(penalty_cov)
  
      penalty_mean = np.asarray(penalty_mean)
      if penalty_mean.shape == ():
          penalty_mean = np.ones(q) * penalty_mean
  
      # compute orthonormal basis of design space
  
      if design is not None:
          X = orth(design)
      else:
          X = None
  
      _iter = 0
      _F = np.inf
  
      while True:
          # Current estimate of mean parameter
  
          iC = npl.inv(C + np.identity(m) / np.exp(32))
  
          # E-step: conditional covariance
  
          if X is not None:
              iCX = np.dot(iC, X)
              Cq = npl.inv(X.T, iCX)
              P = iC - np.dot(iCX, np.dot(Cq, iCX))
          else:
              P = iC
  
          # M-step: ReML estimate of hyperparameters
  
          # Gradient dF/dh (first derivatives)
          # Expected curvature (second derivatives)
  
          U = np.identity(m) - np.dot(P, sigma) / n
  
          for i in range(q):
              PQ[i] = np.dot(P, Q[i])
              dFdh[i] = -(PQ[i] * U).sum() * n / 2
  
              for j in range(i + 1):
                  dFdhh[i, j] = -(PQ[i] * PQ[j]).sum() * n / 2
                  dFdhh[j, i] = dFdhh[i, j]
  
          # Enforce penalties:
  
          dFdh = dFdh - np.dot(penalty_cov, h - penalty_mean)
          dFdhh = dFdhh - penalty_cov
  
          dh = npl.solve(dFdhh, dFdh)
          h -= dh
          C = np.sum([h[i] * Q[i] for i in range(Q.shape[0])], axis=0)
  
          df = (dFdh * dh).sum()
          if np.fabs(df) < 1.0e-01:
              break
  
          _iter += 1
          if _iter >= niter:
              break
  
      return C, h, -dFdhh
  
imports: import numpy
  import numpy.linalg
label_desc: Use image node_id 2 for calling a global function with example usage: reml(sigma, components, design, n, niter, penalty_cov, penalty_mean) and returns: C,  h,  unknown