code: def null_scale(self):
          
          n = self._nsamples
          beta = self.null_beta
          sqrdot = self._yTBy - dot(sum(self._yTBX), beta)
          return sqrdot / n
docs: Optimal s according to the marginal likelihood.
  
          The optimal s is given by ::
  
              s = n⁻¹𝐲ᵀB(𝐲 - X𝜷),
  
          where 𝜷 is optimal.
  
          Returns
          -------
          scale : float
              Optimal scale.