'Estimate sigma. If ML is True, return the ML estimate of sigma,
else return the REML estimate.
If ML, this is (3.6) in Laird, Lange, Stram (see help(Mixed)),
otherwise it corresponds to (3.8).
sigma is the standard deviation of the noise (residual)'

def _compute_sigma(self, ML=False):

    sigmasq = 0.0
    for unit in self.units:
        if ML:
            W = unit.W
        else:
            unit.compute_P(self.Sinv)
            W = unit.P
        t = (unit.r - np.dot(unit.Z, unit.b))
        sigmasq += np.power(t, 2).sum()
        sigmasq += ((self.sigma ** 2) * np.trace((np.identity(unit.n) - ((self.sigma ** 2) * W))))
    self.sigma = np.sqrt((sigmasq / self.N))
