'Approximate covariance of estimates of random effects. Just after
Display (3.10) in Laird, Lange, Stram (see help(Unit)).
D - D\' Z\' P Z D
Notes
In example where the mean of the random coefficient is not zero, this
is not a covariance but a non-centered moment. (proof by example)'

def cov_random(self, D, Sinv=None):

    if (Sinv is not None):
        self.compute_P(Sinv)
    t = np.dot(self.Z, D)
    return (D - np.dot(np.dot(t.T, self.P), t))
