code: def null_beta_covariance(self):
          
          A = sum(i @ j.T for (i, j) in zip(self._XTQDi, self._XTQ))
          return self.null_scale * pinv(A)
docs: Covariance of the optimal 𝜷 according to the marginal likelihood.
  
          Returns
          -------
          beta_covariance : ndarray
              (Xᵀ(s(K + vI))⁻¹X)⁻¹.