code: def bsmahal(a, b, n_boot=200):
      
      n, m = b.shape
      MD = np.zeros((n, n_boot))
      nr = np.arange(n)
      xB = np.random.choice(nr, size=(n_boot, n), replace=True)
  
      # Bootstrap the MD
      for i in np.arange(n_boot):
          s1 = b[xB[i, :], 0]
          s2 = b[xB[i, :], 1]
          X = np.column_stack((s1, s2))
          mu = X.mean(0)
          _, R = np.linalg.qr(X - mu)
          sol = np.linalg.solve(R.T, (a - mu).T)
          MD[:, i] = np.sum(sol**2, 0) * (n - 1)
  
      # Average across all bootstraps
      return MD.mean(1)
docs: Bootstraps Mahalanobis distances for Shepherd's pi correlation.
  
      Parameters
      ----------
      a : ndarray (shape=(n, 2))
          Data
      b : ndarray (shape=(n, 2))
          Data
      n_boot : int
          Number of bootstrap samples to calculate.
  
      Returns
      -------
      m : ndarray (shape=(n,))
          Mahalanobis distance for each row in a, averaged across all the
          bootstrap resamples.