docs: def dirichlet_covariance(alpha):
      r"""Covariance matrix for Dirichlet distribution.
  
      Parameters
      ----------
      alpha : (M, ) ndarray
          Parameters of Dirichlet distribution
  
      Returns
      -------
      cov : (M, M) ndarray
          Covariance matrix
  
      """
      alpha0 = alpha.sum()
      norm = alpha0 ** 2 * (alpha0 + 1.0)
  
      """Non normalized covariance"""
      Z = -alpha[:, np.newaxis] * alpha[np.newaxis, :]
  
      """Correct diagonal"""
      ind = np.diag_indices(Z.shape[0])
      Z[ind] += alpha0 * alpha
  
      """Covariance matrix"""
      cov = Z / norm
  
      return cov