code: def null_scale(self):
          
          np = self._nsamples * self._ntraits
          b = vec(self.null_beta)
          mKiy = b.T @ self._MKiy
          sqrtdot = self._yKiy - mKiy
          scale = sqrtdot / np
          return scale
docs: Optimal s according to the marginal likelihood.
  
          The optimal s is given by
  
              s = (n·p)⁻¹𝐲ᵀK⁻¹(𝐲 - 𝐦),
  
          where 𝐦 = (A ⊗ X)vec(𝚩) and 𝚩 is optimal.
  
          Returns
          -------
          scale : float
              Optimal scale.