code: def est_update_covariance(self, modelparams):
          
          if self._diagonal:
              cov = (self._fixed_scale ** 2 if self._has_fixed_covariance \
                  else np.mean(modelparams[:, self._srw_idxs] ** 2, axis=0))
              cov = np.diag(cov)
          else:
              if self._has_fixed_covariance:
                  cov = np.dot(self._fixed_chol, self._fixed_chol.T)
              else:
                  chol = np.zeros((modelparams.shape[0], self._n_rw, self._n_rw))
                  chol[(np.s_[:],) + self._srw_tri_idxs] = modelparams[:, self._srw_idxs]
                  cov = np.mean(np.einsum('ijk,ilk->ijl', chol, chol), axis=0)
          return cov
docs: Returns the covariance of the gaussian noise process for one 
          unit step. In the case where the covariance is being learned,
          the expected covariance matrix is returned.
          
          :param modelparams: Shape `(n_models, n_modelparams)` shape array
          of model parameters.