Rewriting a given multivariate PDF so that it has the form of multivariate normal

Question

Say, we are given a probability distribution function that we suspect is a multivariate normal. In Mathematica, how can we rewrite such PDF to get its covariance matrix and vector of means?

For instance, we have the following PDF:

$$p(x,y,z) = \frac{\exp \left(\frac{1}{2} \left(\frac{\left(y-x r_{\text{xy}}\right){}^2}{r_{\text{xy}}^2-1}+\frac{z \left(z-x r_{\text{xz}}\right)}{r_{\text{xz}}^2-1}+\frac{x \left(x-z r_{\text{xz}}\right)}{r_{\text{xz}}^2-1}\right)\right)}{2 \sqrt{2} \pi ^{3/2} \sqrt{\left(r_{\text{xy}}^2-1\right) \left(r_{\text{xz}}^2-1\right)}}$$

What is its covariance matrix and the vector of means? How to find them for any PDF that we suspect is a normal distribution?

Answer 1

Here is a partial answer. This answer assumes that there is a polynomial of $x$, $y$, and $z$ in Exp[] and that it has the proper terms such as $x^2$, $y^2$ and $z^2$ and no terms such as $x^3$, $x^2 y$ or Sin[x]. So actually, this partial answer assumes that you already have a multivariate normal pdf and you just want to know the associated mean vector and covariance matrix.

First generate a CoefficientList for the exponential term in the pdf of a general trivariate normal distribution:

o = {x, y, z};
μ = {μx, μy, μz};
Σ = {{σxx, σxy, σxz}, {σxy, σyy, σyz}, {σxz, σyz, σzz}};
n = Expand[-(o - μ).Inverse[Σ].(o - μ)/2];
nc = FullSimplify[CoefficientList[n, o]];

Now do the same thing for your customized distribution:

m = Expand[(1/2) ((y - x rxy)^2/(rxy^2 - 1) + 
      z (z - x rxz)/(rxz^2 - 1) + x (x - z rxz)/(rxz^2 - 1))];
mc = FullSimplify[CoefficientList[m, o]];

Generate a set of equations that matches terms and solve:

eq = Flatten[Table[nc[[i, j, k]] == mc[[i, j, k]], {i, 3}, {j, 3}, {k, 3}]];
sol = Solve[eq, Variables[Flatten[{μ, Σ}]]][[1]]
(* {μx -> 0, μy -> 0, μz -> 0, σxx -> 1, σxy -> rxy, σxz -> rxz,
  σyy -> 1, σyz -> rxy rxz, σzz -> 1} *)

Answer 2

Here's my approach (for what it's worth). Using the definitions of o and m in @Jim Baldwin's answer:

mu = o /. Solve[Thread[D[m, {o}] == 0], o][[1]]    
(* {0, 0, 0} *)

Sigma = Inverse[-D[m, {o, 2}]] // Simplify
(* {{1, rxy, rxz}, {rxy, 1, rxy rxz}, {rxz, rxy rxz, 1}} *)