Multiplying two multivariate Gaussians

[Pi-Bond]

Homework Statement

I am trying to find the resultant Gaussian distribution when two multivariate Gaussians are multiplied together - i.e. find the resultant Fisher matrix and mean.

Homework Equations

Let the two distributions be

P_1(x) = \frac{|A|^{0.5}}{(2\pi)^\frac{n}{2}} exp (-0.5 (x-a)^T A (x-a))
P_2(x) = \frac{|B|^{0.5}}{(2\pi)^\frac{n}{2}} exp (-0.5 (x-b)^T B (x-b))

where A,B are the n-by-n Fisher matrices and a,b are n dimensional mean vectors of the distributions.

The Attempt at a Solution

So I want to find a distribution

P(x) = P_1(x)P_2(x) = P_{0} exp (-0.5 (x-c)^T C (x-c))

where C and c are expressed in terms of A,B,a and b. I've been trying to manipulate the exponents for some time now, but I can't make any progress. Any help would be appreciated.

Thanks.

 

[Ray Vickson]

Homework Statement

I am trying to find the resultant Gaussian distribution when two multivariate Gaussians are multiplied together - i.e. find the resultant Fisher matrix and mean.

Homework Equations

Let the two distributions be

P_1(x) = \frac{|A|^{0.5}}{(2\pi)^\frac{n}{2}} exp (-0.5 (x-a)^T A (x-a))
P_2(x) = \frac{|B|^{0.5}}{(2\pi)^\frac{n}{2}} exp (-0.5 (x-b)^T B (x-b))

where A,B are the n-by-n Fisher matrices and a,b are n dimensional mean vectors of the distributions.

The Attempt at a Solution

So I want to find a distribution

P(x) = P_1(x)P_2(x) = P_{0} exp (-0.5 (x-c)^T C (x-c))

where C and c are expressed in terms of A,B,a and b. I've been trying to manipulate the exponents for some time now, but I can't make any progress. Any help would be appreciated.

Thanks.

Have you forgotten that exp(U)*exp(V) = exp(U+V)?

 

[Pi-Bond]

No, but I can't seem to manipulate the multiplied exponential into the from I want. Just considering the exponent of the multiplied distribution: (-0.5 removed)

(x-a)^TA (x-a) + (x-b)^T B (x-b)
x^TAx - x^TAa -a^TAx +a^TAa + x^TBx -x^TBb -b^tBx +b^TBb

Here a^TAa and b^TBb are just constant so I can absorb them into the constant P_0. I'm not sure what strategy to use on the remaining expression though:

x^TAx - x^TAa -a^TAx + x^TBx -x^tBb -b^tBx

 

[Ray Vickson]

No, but I can't seem to manipulate the multiplied exponential into the from I want. Just considering the exponent of the multiplied distribution: (-0.5 removed)

(x-a)^TA (x-a) + (x-b)^T B (x-b)
x^TAx - x^TAa -a^TAx +a^TAa + x^TBx -b^tBb -b^tBx +b^TBb

Here a^TAa and b^TBb are just constant so I can absorb them into the constant P_0. I'm not sure what strategy to use on the remaining expression though:

x^TAx - x^TAa -a^TAx + x^TBx -b^tBb -b^tBx

Sometimes it is easier to see what is happening by writing things out in detail:
x^T A x = \sum_{i=1}^n \sum_{j=1}^n a_{i,j}\, x_i x_j,
etc.

 

[Pi-Bond]

Ok, so if I apply your suggestion to

x^TAx - x^TAa -a^TAx + x^TBx -x^tBb -b^tBx

I think I should get:

\displaystyle\sum_{ij}^{n} A_{ij}x_{i}x_{j} - A_{ij}x_i a_j -A_{ij}x_i a_j + B_{ij}x_i x_j -B_{ij} x_i b_j -B_{ij} x_i b_j
\displaystyle\sum_{ij}^{n} A_{ij}x_{i}x_{j} - 2A_{ij}x_i a_j + B_{ij}x_i x_j -2B_{ij} x_i b_j

Is this correct? How should I proceed from here? I think I probably need to factor the expression to something like \displaystyle\sum_{ij}^{n} (x_i x_j - ...)(A_{ij}+B_{ij})

 

[Ray Vickson]

Ok, so if I apply your suggestion to

x^TAx - x^TAa -a^TAx + x^TBx -x^tBb -b^tBx

I think I should get:

\displaystyle\sum_{ij}^{n} A_{ij}x_{i}x_{j} - A_{ij}x_i a_j -A_{ij}x_i a_j + B_{ij}x_i x_j -B_{ij} x_i b_j -B_{ij} x_i b_j
\displaystyle\sum_{ij}^{n} A_{ij}x_{i}x_{j} - 2A_{ij}x_i a_j + B_{ij}x_i x_j -2B_{ij} x_i b_j

Is this correct? How should I proceed from here? I think I probably need to factor the expression to something like \displaystyle\sum_{ij}^{n} (x_i x_j - ...)(A_{ij}+B_{ij})

Using the notation <x,Ax> instead of x^TAx (where <u,v> = sum u_iv_i is the inner product) you want to represent
\langle x-a,A(x-a)\rangle+\langle x-b,B(x-b) \rangle
in the form
\langle x-c,C(x-c) \rangle + \,r, where r is a constant. You have already determined that C = A+B, so now you need to determine c.

 

[Pi-Bond]

Well, I don't really know C=A+B from my own calculation. I know it is the answer, but I want to derive it.

&lt;x-a,A(x-a)&gt;+&lt;x-b,B(x-b)&gt;
&lt;x,Ax&gt;-&lt;x,Aa&gt;-&lt;a,Ax&gt;+&lt;a,Aa&gt;+&lt;x,Bx&gt;-&lt;x,Bb&gt;-&lt;b,Bx&gt;+&lt;b,Bb&gt;
&lt;x,Ax&gt;-&lt;x,Aa&gt;-&lt;Ax,a&gt;+&lt;x,Bx&gt;-&lt;x,Bb&gt;-&lt;Bx,b&gt;+const
&lt;x,(A+B)x&gt;-&lt;x,Aa+Bb&gt;-&lt;Ax,a&gt;-&lt;Bx,b&gt;+const

I'm thinking something needs to be added and subtracted here to continue, is that correct?

 

[D H]

I am trying to find the resultant Gaussian distribution when two multivariate Gaussians are multiplied together - i.e. find the resultant Fisher matrix and mean.

Except for trivial cases, the product of two random variables each of which is gaussian is not gaussian. Those trivial cases are where one or both of the random variables has zero variance.

Later on it appears you are looking at the sum of two gaussian RVs. This sum is gaussian if the two random variables are independent.

So which is it, product or sum?

 

[Pi-Bond]

I am trying to show that the product of two multivariate gaussians is also a multivariate gaussian (with another Fisher matrix and mean vector).

The sum which the past few posts show is the exponential part of the multiplied function (see OP). I'm not sure why you are saying the product of two gaussians is not a gaussian. In the univariate case it is true. For example see

https://ccrma.stanford.edu/~jos/sasp/Product_Two_Gaussian_PDFs.html

And I know this can be generalized to the multivariate case - I'm working from a question which says it is. I'm just having troubles trying to prove it!

 

[Pi-Bond]

I think I got it! If I expand
&lt;x-c,C(x-c)&gt;= &lt;x,Cx&gt;-&lt;c,Cx&gt;-&lt;x,Cc&gt;+&lt;c,Cc&gt;
And compare with the previous expansion, I find:
C=A+B
c=(A+B)^{-1} Aa+(A+B)^{-1}Bb

Thanks a lot for the help, Ray Vickson, it is much appreciated.

 

[D H]

For example see

https://ccrma.stanford.edu/~jos/sasp/Product_Two_Gaussian_PDFs.html

And I know this can be generalized to the multivariate case - I'm working from a question which says it is. I'm just having troubles trying to prove it!

Everything on the internet is true!

Except when it isn't. Let's look at this *bad* math from a number of perspectives, analytically, finding cases where this fails, and Monte Carlo simulation.

1. Analytically.
By definition Cov[x,y] = E[(x-E(x))*(y-E(y)]. Expanding this, one gets Cov[x,y] = E[x*y] -E[x]*E[y]. In other words, E[x*y] = Cov[x,y] + E[x]*E[y]. In words, the expected value of the product of two one dimensional random variables is the sum of the covariance and the product of the expected values. The formula that you found for the mean (and yes, it's all over the internet) doesn't look anything like this.

2. Find obvious cases where the formulae are wrong.
Case 1: σ_y is zero (in other words, y is constant). A constant times a gaussian is a guassian, but the mean and variance are not anywhere close to the values given by those formulae.
Case 2: Let Y=X (correlation=1). Now X*Y is always non-negative - so it can't be gaussian.

3. Monte Carlo simulation.
Here's a simple python script.

import math
import random

N = 100000

mu_x  = 20
sig_x = 20

mu_y  = 2 
sig_y = 10

x = [random.gauss(mu_x,sig_x) for ii in range(N)]
y = [random.gauss(mu_y,sig_y) for ii in range(N)]

def prod (X,Y) : return X*Y 
z = map (prod, x, y)def report (name, mu, sig, X) :
   N = len(X)
   mean = sum(X) / N 

   xsq = 0.0 
   for val in X : 
      xsq = xsq + (val - mean)**2
   stddev = math.sqrt(xsq/(N-1))

   print "\n" + name + ":" 
   print "expected mean, std_dev " + str(mu) + ", " + str(sig)
   print "observed mean, std_dev " + str(mean) + ", " + str(stddev)report ("x", mu_x, sig_x, x)
report ("y", mu_y, sig_y, y)

mu_z  = (mu_x*sig_y**2 + mu_y*sig_x**2)/(sig_x**2 + sig_y**2)
sig_z = math.sqrt(((sig_x**2) * (sig_y**2)) / (sig_x**2 + sig_y**2))
report ("z", mu_z, sig_z, z)

The product of two gaussians is not a gaussian except in the trivial case that one or both of them has zero variance.

 

[Ray Vickson]

Everything on the internet is true!

Except when it isn't. Let's look at this *bad* math from a number of perspectives, analytically, finding cases where this fails, and Monte Carlo simulation.

1. Analytically.
By definition Cov[x,y] = E[(x-E(x))*(y-E(y)]. Expanding this, one gets Cov[x,y] = E[x*y] -E[x]*E[y]. In other words, E[x*y] = Cov[x,y] + E[x]*E[y]. In words, the expected value of the product of two one dimensional random variables is the sum of the covariance and the product of the expected values. The formula that you found for the mean (and yes, it's all over the internet) doesn't look anything like this.

2. Find obvious cases where the formulae are wrong.
Case 1: σ_y is zero (in other words, y is constant). A constant times a gaussian is a guassian, but the mean and variance are not anywhere close to the values given by those formulae.
Case 2: Let Y=X (correlation=1). Now X*Y is always non-negative - so it can't be gaussian.

3. Monte Carlo simulation.
Here's a simple python script.

import math
import random

N = 100000

mu_x  = 20
sig_x = 20

mu_y  = 2 
sig_y = 10

x = [random.gauss(mu_x,sig_x) for ii in range(N)]
y = [random.gauss(mu_y,sig_y) for ii in range(N)]

def prod (X,Y) : return X*Y 
z = map (prod, x, y)


def report (name, mu, sig, X) :
   N = len(X)
   mean = sum(X) / N 

   xsq = 0.0 
   for val in X : 
      xsq = xsq + (val - mean)**2
   stddev = math.sqrt(xsq/(N-1))

   print "\n" + name + ":" 
   print "expected mean, std_dev " + str(mu) + ", " + str(sig)
   print "observed mean, std_dev " + str(mean) + ", " + str(stddev)


report ("x", mu_x, sig_x, x)
report ("y", mu_y, sig_y, y)

mu_z  = (mu_x*sig_y**2 + mu_y*sig_x**2)/(sig_x**2 + sig_y**2)
sig_z = math.sqrt(((sig_x**2) * (sig_y**2)) / (sig_x**2 + sig_y**2))
report ("z", mu_z, sig_z, z)

The product of two gaussians is not a gaussian except in the trivial case that one or both of them has zero variance.

He is not multiplying the random variables. He is multiplying two Gaussian functions. I did not check whether this is a sensible thing to do; I just took him at his word.

 

[Pi-Bond]

It would be a tad odd if Stanford was concocting falsehoods! I did post with the intent of multiplying two gaussian functions rather than two gaussian distributed variables. The context is in Bayes Theorem; one gaussian is the prior, the other the likelihood. What are your thoughts on the matter?

Also, I am actually approaching the matter from a Physics point of view. I won't be surprised if the mathematical rigour is not upto mark in our calculations.

 

[D H]

That's different. I thought you were taking about the product of two gaussian RVs. The product of two gaussian functions is indeed another gaussian (but note: it's no longer normalized).

 

[Pi-Bond]

Alright, thanks for the clarification.