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?


2 Answers 2


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} *)
  • $\begingroup$ Nice! This solves my question. @Jim writes that this is a "partial solution", but a more generic solution would be solving a harder more generic problem, stated in this stackexchange post $\endgroup$
    – pms
    Aug 11, 2017 at 17:20

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}} *)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.