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

Learn moreOpportunities for recent engineering grads.

Apply Today
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

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

## 11 Comments

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?

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

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.

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.

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 ?

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

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!

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.

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.

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.

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

Nice answer! Thank a lot!