code: def estimate_pooled_covariance(resid, ARtarget=[0.3], mask=None):
      """
      Use SPM's REML implementation to estimate a pooled covariance matrix.
  
      Thresholds an F statistic at a marginal pvalue to estimate
      covariance matrix.
  
      """
      resid
      n = resid[:].shape[0]
      components = correlation.ARcomponents(ARtarget, n)
  
      raw_sigma = 0
      nvox = 0
      for i in range(resid.shape[1]):
          d = np.asarray(resid[:, i])
          d.shape = (d.shape[0], np.prod(d.shape[1:]))
          keep = np.asarray(mask[i])
          keep.shape = np.prod(keep.shape)
          d = d.compress(keep, axis=1)
          raw_sigma += np.dot(d, d.T)
          nvox += d.shape[1]
      raw_sigma /= nvox
      C, h, _ = reml.reml(raw_sigma, components, n=nvox)
      return C
  
file_path: nipy/nipy/modalities/fmri/spm/model.py
imports: import numpy
  import numpy.linalg
  from scipy.stats import f
  from nipy.algorithms.statistics.models.regression import GLSModel, OLSModel
  from nipy.core.api import Image
  from nipy.modalities.fmri.fmristat import model
  from nipy.modalities.fmri.fmristat.model import OLS
  from .None import correlation, reml
label_desc: Use image node_id 2 for calling a global function with example usage: estimate_pooled_covariance(resid, ARtarget, mask) and returns: C