code: def posterior(self, x, Sigma=None):
          
          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)
docs: Posterior for model: X1, ..., Xn ~ N(theta, Sigma).
          
          Parameters
          ----------
          x: (n, d) ndarray
              data 
          Sigma: (d, d) ndarray 
              (fixed) covariance matrix in the model