How to calculate the product of Gaussian variables in Matlab

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?

• Is your normal equation right? I would have thought: normal = @(x, mu, Sig) (1/(Sig*sqrt(2*pi))) * exp((-1/2)*((x-mu)/Sig)^2);? Not that it affects your charts – Dan Jun 25 '14 at 9:13
• You're right, it looked a bit confusing. I changed Sig to sig2 to make clear I use sigma^2 in my definition. – Edgar Jun 25 '14 at 11:37

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.

• Thanks for your answer, Tom! I tried it, but it still doesn't match, see my updated question. If you have further ideas, please let me know! – Edgar Nov 25 '14 at 15:29
• I left out an absolute value in the Jacobian. It should be correct now. – Tom Minka Nov 26 '14 at 10:43
• Thanks, Tom! That did the trick! Now the numeric integral matches the sampling solution pretty well. (Anyway, the Bessel-version still doesn't match.) – Edgar Nov 26 '14 at 11:01

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.

• Good point! I flagged my question and suggested a move to Cross Validated. – Edgar Jul 1 '14 at 10:17

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')


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.

• Welcome to our site. When you believe that all answers to a question are "completely wrong" (which does happen some of the time--roughly in about 0.01% of cases), it is wise to take that to mean you are interpreting the question differently than other respondents. To avoid being misunderstood (and severely downvoted) it would be a great idea to clarify your explanation of just how you are interpreting the question. I'm having a hard time thinking of any other way to understand "product of two iid Gaussian distributed variables," but I'm ready to learn. – whuber Sep 20 '18 at 21:30
• Here is the correct Matlab codes and it works: – Sahak Khacheryan Sep 27 '18 at 22:49