How to calculate the product of Gaussian variables in Matlab

Question

I want to calculate the distribution of a product of two i.i.d. Gaussian distributed variables a and b. In principle, this should be possible by defining a new variable x with a dirac delta distribution

To get the distribution over x (the product of a and b), a and b have to be marginalized out:

Since the dirac delta ties the variables a and b together via the variable x=ab, we can skip one integral and get

So far so good.

If I now want to calculate this in Matlab, I get a huge difference between the solution with integral and a simple sampling solution as approximation.

Here is my MWE:

mu_a = 2;
sig2_a = 1;

mu_b = 0;
sig2_b = 1;

res = 250;
xmin = -10;
xmax = 10;
x_range = linspace(xmin,xmax,res);

normal = @(x, mu, sig2) 1/sqrt(2*pi*sig2)*...
    exp((-1/2)*((x-mu)^2/sig2));

normal_a = @(a) arrayfun(normal,a,mu_a*ones(size(a)),sig2_a*ones(size(a)));
normal_b = @(b) arrayfun(normal,b,mu_b*ones(size(b)),sig2_b*ones(size(b)));

product_pdf1 = zeros(size(x_range));
product_pdf2 = zeros(size(x_range));
for i = 1:length(x_range)    
    x = x_range(i);

    normal_ab = @(a) normal_a(a).*normal_b(x./a);
    product_pdf1(i) = integral(normal_ab,-Inf,Inf);

    normal_ab = @(b) normal_a(x./b).*normal_b(b);
    product_pdf2(i) = integral(normal_ab,-Inf,Inf);
end

subplot(2,1,1,'replace'); hold all
plot(x_range, normal_a(x_range))
plot(x_range, normal_b(x_range))
plot(x_range, product_pdf1, '--')
plot(x_range, product_pdf2, '--')
axis tight
drawnow

% the sampling solution
N = 50000;
X_sequence = zeros(1,N);
for i=1:N
    a = mu_a + sqrt(sig2_a)*randn();
    b = mu_b + sqrt(sig2_b)*randn();
    x = a*b;
    X_sequence(i) = x;
end

subplot(2,1,2,'replace'); hold all
hist(X_sequence,250)
xlim([xmin xmax])

Note that I calculated the integral once for a and once for b, the output is completely different.

Depending on the parameters, one, the other or both of the calculations are clearly wrong. In my example, the sampling solution is much peakier than the better looking integral attempt. The other one is completely wrong.

I now wonder: Is my math wrong? Is my code wrong? Or is it a limitation of Matlab's integral function (e.g. having trouble integrating over discontinuities)?

UPDATE:

Tom spotted a mistake in my derivations. I updated my MWE accordingly and normalized the histogram for better comparison:

mu_a = 2;
sig2_a = 1;

mu_b = 0;
sig2_b = 1;

res = 250;
xmin = -10;
xmax = 10;
x_range = linspace(xmin,xmax,res);

normal = @(x, mu, sig2) 1/sqrt(2*pi*sig2)*...
    exp((-1/2)*((x-mu)^2/sig2));

normal_a = @(a) arrayfun(normal,a,mu_a*ones(size(a)),sig2_a*ones(size(a)));
normal_b = @(b) arrayfun(normal,b,mu_b*ones(size(b)),sig2_b*ones(size(b)));

product_pdf1 = zeros(size(x_range));
product_pdf2 = zeros(size(x_range));
for i = 1:length(x_range)    
    x = x_range(i);

    normal_ab = @(a) (1./a).*normal_a(a).*normal_b(x./a);
    product_pdf1(i) = integral(normal_ab,-Inf,Inf);

    normal_ab = @(b) (1./b).*normal_a(x./b).*normal_b(b);
    product_pdf2(i) = integral(normal_ab,-Inf,Inf);
end

subplot(2,1,1,'replace'); hold all
plot(x_range, normal_a(x_range))
plot(x_range, normal_b(x_range))
plot(x_range, product_pdf1, '--')
plot(x_range, product_pdf2, '--')
axis tight
drawnow

% the sampling solution
N = 50000;
X_sequence = zeros(1,N);
for i=1:N
    a = mu_a + sqrt(sig2_a)*randn();
    b = mu_b + sqrt(sig2_b)*randn();
    x = a*b;
    X_sequence(i) = x;
end

subplot(2,1,2,'replace'); hold all
hres = 250; % histogram resolution
[B, X] = hist(X_sequence,hres);
dx = diff(X(1:2));
plot(x_range, product_pdf1,'--')
plot(X, B/sum(B*dx),'-.');
plot(x_range, besselk(0,abs(x_range))/pi)
xlim([xmin xmax])

But, strangely, the solutions still don't match: One of the integral solutions is now completely off, but the other (lower plot, dot-dashed red) is now actually closer to the sampling solution (dashed green). The BesselK doesn't match, it is much peakier.

Further suggestions?

Answer 1

You made a mistake in the derivation. Let's start from this integral: $$ p(x) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \delta(ab-x) p(a) p(b) da ~ db $$ To remove the $\delta$, you need to make a change of variable from $b$ to $c=ab$, for which $dc = a ~ db$. The integral becomes: $$ p(x) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \delta(c-x) p(a) p(c/a) \frac{1}{|a|} da ~ dc \\ = \int_{-\infty}^{\infty} p(a) p(x/a) \frac{1}{|a|} da $$ When $p(a)$ is standard normal, $p(x) = \frac{1}{\pi} {\rm BesselK}(0,|x|)$. You can compute this in Matlab with the code besselk(0,abs(x))/pi.

Answer 2

Your MATLAB code looks fine. It also agrees with my own implementation of your equations:

mu2 = 2;
[a, x] = ndgrid(linspace(-10, 10, 250));
clf; hold all
hist(randn(1,1e6) .* (randn(1, 1e6) + mu2), x(1,:));
y = sum(exp(-0.5.*(a-mu2).^2).*exp(-0.5.*(x./a).^2))';
y = y * (1e6 / sum(y));
plot(x(1,:)', y, 'r-');
y = sum(exp(-0.5.*a.^2).*exp(-0.5.*((x./a)-mu2).^2))';
y = y * (1e6 / sum(y));
plot(x(1,:)', y, 'g-');
xlim(x(1,[1 end]));

The fact that our two different approaches to computing the integrals agree, and that the two integrals disagree both with each other and with the empirical distribution, suggests to me that the problem lies with your mathematical derivation. However, I couldn't spot it, but that question isn't appropriate for this exchange anyway. I suggest you post it on the Mathematics or Cross Validated exchanges.

Answer 3

I think your derivation is messed up. Take a look at this Wiki.

Here's MATLAB implementation:

 f1=@(x,z)normpdf(x).*normpdf(z./x)./abs(x)
 f=@(z)integral(@(x)f1(x,z),-10,10)
 >> f(1)
 ans =
 0.1340

Here's how you could simulate it:

 samp=prod(randn(2,100000));
 hist(samp,100);

Use distribution fitting tool to fit nonparametric distribution to compare to theoretical distribution:

 dfittool(samp')

Answer 4

The sampling solution above in Matlab codes is completely wrong. pdf function of Z=X*Y means that X and Y values could be anything as long as their product is equal to Z!!!!

For correct and computationally optimal answer please see paper by: http://citeseerx.ist.psu.edu/viewdoc/downloaddoi=10.1.1.15.7546&rep=rep1&type=pdf Andrew G. Glen, Lawrence M. Leemis, John H., Drew, "Computing the Distribution of the Product of Two Continuous Random Variables," Computational Statistics & Data Analysis, 2004. In this paper the lower limits of PDF functions a & b are larger than 0. I hope it helps.