code: def _mc_error(x, batches=5, circular=False):
      
      if x.ndim > 1:
  
          dims = np.shape(x)
          trace = np.transpose([t.ravel() for t in x])
  
          return np.reshape([_mc_error(t, batches) for t in trace], dims[1:])
  
      else:
          if batches == 1:
              if circular:
                  std = st.circstd(x, high=np.pi, low=-np.pi)
              else:
                  std = np.std(x)
              return std / np.sqrt(len(x))
  
          batched_traces = np.resize(x, (batches, int(len(x) / batches)))
  
          if circular:
              means = st.circmean(batched_traces, high=np.pi, low=-np.pi, axis=1)
              std = st.circstd(means, high=np.pi, low=-np.pi)
          else:
              means = np.mean(batched_traces, 1)
              std = np.std(means)
  
          return std / np.sqrt(batches)
docs: Calculate the simulation standard error, accounting for non-independent samples.
  
      The trace is divided into batches, and the standard deviation of the batch
      means is calculated.
  
      Parameters
      ----------
      x : Numpy array
          An array containing MCMC samples
      batches : integer
          Number of batches
      circular : bool
          Whether to compute the error taking into account `x` is a circular variable
          (in the range [-np.pi, np.pi]) or not. Defaults to False (i.e non-circular variables).
  
      Returns
      -------
      mc_error : float
          Simulation standard error