code: def _calc_ARMA_noise(
      volume,
      mask,
      auto_reg_order=1,
      ma_order=1,
      sample_num=100,
  ):
      """Calculate the the ARMA noise of a volume
      This calculates the autoregressive and moving average noise of the volume
      over time by sampling brain voxels and averaging them.
  
      Parameters
      ----------
  
      volume : 4d array or 1d array, float
          Take a volume time series to extract the middle slice from the
          middle TR. Can also accept a one dimensional time course (mask input
          is then ignored).
  
      mask : 3d array, binary
          A binary mask the same size as the volume
  
      auto_reg_order : int
          What order of the autoregression do you want to estimate
  
      sample_num : int
          How many voxels would you like to sample to calculate the AR values.
          The AR distribution of real data is approximately exponential maxing
          at 1. From analyses across a number of participants, to get less
          than 3% standard deviation of error from the true mean it is
          necessary to sample at least 100 voxels.
  
      Returns
      -------
  
      auto_reg_rho : list of floats
          Rho of a specific order for the autoregression noise in the data
  
      na_rho : list of floats
          Moving average of a specific order for the data
  
      """
  
      # Pull out the non masked voxels
      if len(volume.shape) > 1:
          brain_timecourse = volume[mask > 0]
      else:
          # If a 1 dimensional input is supplied then reshape it to make the
          # timecourse
          brain_timecourse = volume.reshape(1, len(volume))
  
      # Identify some brain voxels to assess
      voxel_idxs = list(range(brain_timecourse.shape[0]))
      np.random.shuffle(voxel_idxs)
  
      # If there are more samples than voxels, take all of the voxels
      if len(voxel_idxs) < sample_num:
          sample_num = len(voxel_idxs)
  
      auto_reg_rho_all = np.zeros((sample_num, auto_reg_order))
      ma_all = np.zeros((sample_num, ma_order))
      for voxel_counter in range(sample_num):
          # Get the timecourse and demean it
          timecourse = brain_timecourse[voxel_idxs[voxel_counter], :]
          demeaned_timecourse = timecourse - timecourse.mean()
  
          # Pull out the ARMA values (depends on order)
          try:
              model = ARIMA(
                  demeaned_timecourse,
                  order=[auto_reg_order, 0, ma_order],
              )
              model_fit = model.fit()
              params = model_fit.params
          except (ValueError, LinAlgError):
              params = np.ones(auto_reg_order + ma_order + 1) * np.nan
  
          auto_reg_rho_all[voxel_counter, :] = params[
              1 : auto_reg_order + 1
          ]
          ma_all[voxel_counter, :] = params[auto_reg_order + 1 : -1]
  
      # Average all of the values and then convert them to a list
      auto_reg_rho = np.nanmean(auto_reg_rho_all, 0).tolist()
      ma_rho = np.nanmean(ma_all, 0).tolist()
  
      # Return the coefficients
      return auto_reg_rho, ma_rho
  
imports: import logging
  from itertools import product
  from statsmodels.tsa.arima.model import ARIMA
  import math
  import numpy
  from numpy.linalg import LinAlgError
  from pkg_resources import resource_stream
  from scipy import stats
  from scipy import signal
  import scipy.ndimage
  import copy
  from scipy import optimize
label_desc: Use image node_id 13 for calling a global function with example usage: _calc_ARMA_noise(volume, mask, auto_reg_order, ma_order, sample_num) and returns: auto_reg_rho,  ma_rho