1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Multiplying two multivariate Gaussians

  1. Nov 26, 2012 #1
    1. The problem statement, all variables and given/known data
    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.

    2. Relevant equations
    Let the two distributions be

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

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

    3. The attempt at a solution
    So I want to find a distribution

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

    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.

  2. jcsd
  3. Nov 26, 2012 #2

    Ray Vickson

    User Avatar
    Science Advisor
    Homework Helper

    Have you forgotten that exp(U)*exp(V) = exp(U+V)?
  4. Nov 26, 2012 #3
    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)

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

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

    [tex]x^TAx - x^TAa -a^TAx + x^TBx -x^tBb -b^tBx[/tex]
    Last edited: Nov 26, 2012
  5. Nov 26, 2012 #4

    Ray Vickson

    User Avatar
    Science Advisor
    Homework Helper

    Sometimes it is easier to see what is happening by writing things out in detail:
    [tex] x^T A x = \sum_{i=1}^n \sum_{j=1}^n a_{i,j}\, x_i x_j, [/tex]
  6. Nov 26, 2012 #5
    Ok, so if I apply your suggestion to

    [tex]x^TAx - x^TAa -a^TAx + x^TBx -x^tBb -b^tBx[/tex]

    I think I should get:

    [tex]\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 [/tex]
    [tex]\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[/tex]

    Is this correct? How should I proceed from here? I think I probably need to factor the expression to something like [tex]\displaystyle\sum_{ij}^{n} (x_i x_j - ...)(A_{ij}+B_{ij})[/tex]
  7. Nov 27, 2012 #6

    Ray Vickson

    User Avatar
    Science Advisor
    Homework Helper

    Using the notation <x,Ax> instead of xTAx (where <u,v> = sum uivi is the inner product) you want to represent
    [tex] \langle x-a,A(x-a)\rangle+\langle x-b,B(x-b) \rangle [/tex]
    in the form
    [tex] \langle x-c,C(x-c) \rangle + \,r,[/tex] where r is a constant. You have already determined that C = A+B, so now you need to determine c.
  8. Nov 27, 2012 #7
    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.


    I'm thinking something needs to be added and subtracted here to continue, is that correct?
  9. Nov 27, 2012 #8

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    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?
  10. Nov 27, 2012 #9
    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


    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!
  11. Nov 27, 2012 #10
    I think I got it! If I expand
    [tex]<x-c,C(x-c)>= <x,Cx>-<c,Cx>-<x,Cc>+<c,Cc>[/tex]
    And compare with the previous expansion, I find:
    [tex]c=(A+B)^{-1} Aa+(A+B)^{-1}Bb[/tex]

    Thanks a lot for the help, Ray Vickson, it is much appreciated.
  12. Nov 27, 2012 #11

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    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.
    Code (Text):
    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.
  13. Nov 27, 2012 #12

    Ray Vickson

    User Avatar
    Science Advisor
    Homework Helper

    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.
  14. Nov 27, 2012 #13
    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.
  15. Nov 28, 2012 #14

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    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).
  16. Nov 28, 2012 #15
    Alright, thanks for the clarification.
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Similar Discussions: Multiplying two multivariate Gaussians