Asked by may
on 16 Sep 2013

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.

Answer by Sean de Wolski
on 16 Sep 2013

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

Answer by Walter Roberson
on 16 Sep 2013

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
on 16 Sep 2013

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi test

Learn moreOpportunities for recent engineering grads.

Apply Today
## 11 Comments

## Roger Stafford (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/87362#comment_169286

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 (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/87362#comment_169289

@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 (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/87362#comment_169336

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:

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 (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/87362#comment_169342

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.

## Walter Roberson (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/87362#comment_169344

What error does it give?

What value would you be expecting at x - 0 ?

## may (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/87362#comment_169347

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

## may (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/87362#comment_169350

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

## Roger Stafford (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/87362#comment_169355

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 (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/87362#comment_169528

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 (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/87362#comment_169530

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

## may (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/87362#comment_169562

Nice answer! Thank a lot!