def posterior(self, x, Sigma=None):
        """Posterior for model: X1, ..., Xn ~ N(theta, Sigma).
        
        Parameters
        ----------
        x: (n, d) ndarray
            data 
        Sigma: (d, d) ndarray 
            (fixed) covariance matrix in the model
        """
        n = x.shape[0]
        Sigma = np.eye(self.dim) if Sigma is None else Sigma
        Siginv = inv(Sigma)
        Qpost = inv(self.cov) + n * Siginv
        Sigpost = inv(Qpost)
        mupost = (np.matmul(Siginv, self.mean) + 
                  np.matmu(Siginv, np.sum(x, axis=0)))
        return MvNormal(loc=mupost, cov=Sigpost)