code: def covariance(self):
          r
          from numpy_sugar.linalg import ddot, sum2diag
  
          Q0 = self._QS[0][0]
          S0 = self._QS[1]
          return sum2diag(dot(ddot(Q0, self.v0 * S0), Q0.T), self.v1)
docs: r"""Covariance of the prior.
  
          Returns
          -------
          :class:`numpy.ndarray`
              :math:`v_0 \mathrm K + v_1 \mathrm I`.