@Andrew M You are right, the question was misleading. I've made corresponding edits.

@Andrew M Ideally, yes, but a single solution will also do.

@Dilip Sarwate I understand that, I'm just surprised that it makes such a big difference. I was wondering if it is possible to introduce an uninformative 'dummy' dimension in the one-dimensional density to make Eq.(371) work.

Thank you very much, that helped me a lot! I'm still confused how to map those results to Eq. (371) in the Matrix Cookbook, which says that the product of two Gaussians is proportional to a Gaussian even when the means and covariances are not consistent. The derivation above seems to assume that there is only one $\sigma_x$. In my case, however, the variances of $X$ are different in the one- and two-dimensional density.

Note that the variables in the two densities are $not$ disjunct, so the covariance of $p$ is $not$ block diagonal -- a case which is discussed here

I guess the estimate is the Gaussian with minimal KL-divergence to the true generating distribution. The MSE/bias/variance could give a hint how wrong the Gaussian model assumption is.