'Return log-likelihood, REML by default.'

def logL(self, ML=False):

    logL = 0.0
    for unit in self.units:
        logL += unit.logL(a=self.a, ML=ML)
    if (not ML):
        logL += (np.log(L.det(self.Sinv)) / 2)
    return logL
