How to find double integral in MATLAB

Given mean1, mean2, sigma1,sigma2, and u, I want to find the following integral:

for example: mean1=0, mean2=0, sigma1=0.2, sigma2=1, and u=0.4

F = @(x,y)normpdf(x, mean1, sigma1).*normpdf(y, mean2, sigma2).*dirac(x*y-u);
answer= integral2(F,-Inf,Inf,-Inf,Inf);

but I get the following error

|Error using integralCalc/finalInputChecks (line 511) Input function must return 'double' or 'single' values. Found 'function_handle'.|

I would appreciate if you could help me to fix my code. Thank you.

[Roger Stafford]

You should not be using the 'dirac' function as part of the integrand in a numerical integration function like 'integral2'. The answer obtained by Sean is correct, namely, zero. The 2D integral along the 1D curve, x*y == 0.4, is bound to be zero for otherwise continuous integrands. You did not need to call on 'integral2' to find that out. Are you sure this is the integration you wished to perform?

[may]

@Roger Stafford 0.4 is just an example, the integral i want to find can be found in the following link http://mathworld.wolfram.com/NormalProductDistribution.html I would appreciate if you could help me

[Roger Stafford]

I was mistaken when I said Sean's answer of zero was correct. It is not correct. However, you cannot solve such a double integral involving the dirac delta function using 'integral2' directly for the double integral you saw on the wolfram website. That lies far beyond the capabilities of an ordinary numerical double integration function. As Walter states, you might conceivably obtain good answers using the 'int' of the Symbolic Toolbox, (if you are very lucky.)

There is another conceivable way of solving it which would require only a single numerical integral for each value of u, but I confess it is a bit speculative on my part and depends on one's interpretation of the P(u) in the Wolfram site. If we assume that P(u) is to be a true density distribution with respect to the product parameter u = x*y, then P(u) would be equal to the line integral followed along the two disjoint paths for a constant product u = x*y taken with respect to path arclength using the following integrand:

 f1(x)*f2(y)/sqrt(x^2+y^2)

where x and y are functions of arclength s along the curve and where f1(x) and f2(y) are the two normal distributions in question. (The sqrt(x^2+y^2) divisor is the Jacobian correction needed to make P(u) a true density with respect to the product u.) I won't attempt to tell you how to proceed in carrying out this line integral possibility. It would take some careful thinking. At first thought it might be done using one of matlab's ode solvers.

[may]

Thanks a lot for your explanation and help. I could finally simplify the double integral, now, I should find the following single integral. but when I run it it gives error! I do not know how to fix it. Any help would be appreciated.

   F = @(x)normpdf(x, mean1, sigma1).*normpdf(u./x, mean2, sigma2).*(1./abs(x));
  answer= integral(F,-Inf,Inf);

[Walter Roberson]

What error does it give?

What value would you be expecting at x - 0 ?

[may]

Warning: Infinite or Not-a-Number value encountered. I think as you said the problem is at x=0

[may]

I think it does not work for zero means, but when I use non-zero means it returns some answer!

[Roger Stafford]

Yes, that form is simpler and it looks valid to me. The price you pay, of course, is the computational difficulties in the neighborhood of x equal to zero. It takes on the nature of (almost) zero divided by zero there which could produce NaNs. You may have to pull short of x = 0 a small amount on either side to avoid this. L'Hospital's rule shows that x = 0 is not actually a singularity, but that doesn't mean 'integral' won't have trouble there.

[may]

Thanks for your help. I think another way is to approximate the distribution of X*Y by producing some random numbers from the distribution of X and Y then multiply them all together and find its histogram.

[Walter Roberson]

If you define a true function you can use conditional logic for it.

function r = F(x)
  r = zeros(size(x));
  r(x==0) = whatever it should be;
  nzx = x(x ~=0);
  r(x~=0) = normpdf(nxz, mean1, sigma1).*normpdf(u./nzx, mean2, sigma2).*(1./abs(nzx);
end

[may]

Nice answer! Thank a lot!

Answer 1

[Sean de Wolski]

Running your example I get an error saying the matrix dimensions must agree. Changing the x*y to x.*y inside of dirac allows it to run and returns an answer of zero

F = @(x,y)normpdf(x, mean1, sigma1).*normpdf(y, mean2, sigma2).*dirac(x.*y-u);
answer= integral2(F,-Inf,Inf,-Inf,Inf);

[may]

thank you, but now I get zero for all the inputs i am trying.

Answer 2

[Walter Roberson]

Which dirac are you using? It appears you might be using dirac from the Symbolic toolbox, You should experiment with

class(dirac(5))   %for example

You might need something like double(dirac(x*y-u)) in your code instead of a plan dirac() call.

[Walter Roberson]

You should be using symbolic integration instead of numeric integration. symbolic integration is int() in the Symbolic Toolbox