code: def recursive_cov(self, cov, length, mean, chain, scaling=1, epsilon=0):
          r
          n = length + len(chain)
          k = length
          new_mean = self.recursive_mean(mean, length, chain)
  
          t0 = k * np.outer(mean, mean)
          t1 = np.dot(chain.T, chain)
          t2 = n * np.outer(new_mean, new_mean)
          t3 = epsilon * np.eye(cov.shape[0])
  
          new_cov = (
              k - 1) / (
                  n - 1.) * cov + scaling / (
                      n - 1.) * (
                          t0 + t1 - t2 + t3)
          return new_cov, new_mean
docs: r"""Compute the covariance recursively.
  
          Return the new covariance and the new mean.
  
          .. math::
              C_k & = \frac{1}{k-1} (\sum_{i=1}^k x_i x_i^T - k\bar{x_k}\bar{x_k}^T)
              C_n & = \frac{1}{n-1} (\sum_{i=1}^k x_i x_i^T + \sum_{i=k+1}^n x_i x_i^T - n\bar{x_n}\bar{x_n}^T)
                  & = \frac{1}{n-1} ((k-1)C_k + k\bar{x_k}\bar{x_k}^T + \sum_{i=k+1}^n x_i x_i^T - n\bar{x_n}\bar{x_n}^T)
  
          :Parameters:
              -  cov : matrix
                  Previous covariance matrix.
              -  length : int
                  Length of chain used to compute the previous covariance.
              -  mean : array
                  Previous mean.
              -  chain : array
                  Sample used to update covariance.
              -  scaling : float
                  Scaling parameter
              -  epsilon : float
                  Set to a small value to avoid singular matrices.