I would like to take double integral over standard normal distribution but it does not work.

Hi everyone, I'm trying to take double integral over the pdf to find out the probablity,

However using the code below, it just took the time and memory(over 9GB) without any error message, but it did not give any solution.

What I'm trying: where .

syms z_1 z_2 rho
phi(z_1,z_2,rho) = (exp(-(((z_1)^2)+((z_2)^2)-(2*rho*(z_1)*(z_2)))/(2*(1-rho^2)))/(2*pi*sqrt(1-rho^2)))
 
phi(z_1, z_2, rho) =
 
exp((z_1^2 - 2*rho*z_1*z_2 + z_2^2)/(2*rho^2 - 2))/(2*pi*(1 - rho^2)^(1/2))
 
syms mu_1 mu_2 sigma lambda mu_M
Pr1 = int(int(phi,z_2,[-inf ((mu_M-mu_2)/sigma)]),z_1,[-inf ((lambda-mu_1)/sigma)])

I would appreciate your help.

[Alan Stevens]

There is no symbolic solution when integrating the normal curve (with one or more finite, non-zero limits). You can only do this numerically.

[Bjorn Gustavsson]

@Alan: Yes there is, we get the error-function:

>> syms x xlim real
>> F = int(exp(-x^2),-inf,xlim)
 
F =
 
(pi^(1/2)*(erf(xlim) + 1))/2

@OP: You will find at least some relevant material here: Report about error-function in 2-D

[Alan Stevens]

Ok, so we get the integral in terms of another function which itself is solved numerically.

What I mean is int(exp(-x^2),-inf,xlim) is, by definition, the (complementary) error function between those limits (apart from a normalizing constant), so going from the integral of exp(-x^2) to the error function is mathematically not much more than a tautology.

I guess one can call it progress in MATLAB terms, given that MATLAB has the error function built in, and can evaluate it numerically!

[Bjorn Gustavsson]

Pretty much the same line of reasoning can be made for elliptical integral-functions, Bessel-functions and even the exponential function. The symbolic toolbox can for example differentiate, integrate and produce Taylor-expansions of the error-function, I see no urgent need to stop the symbolic work because the definition of a function becomes to complicated - at least when we have a computer doing the hard work...

[Bruno Luong]

I don't think 2D integral formula of a Gaussian can be expressed as 1D integral and error function.

[Bjorn Gustavsson]

@Bruno, For the uncorrelated case things behave nicely:

>> syms x z xlim zlim real
>> TwoDtrivial = int(int(exp(-x^2-z^2),'x',-inf,xlim),'z',-inf,zlim)
 
TwoDtrivial =
 
(pi*(erf(xlim) + 1)*(erf(zlim) + 1))/4

For the case I sketched below (in its simplest version where the transformation is just a 45-degree rotation) I end up with:

>> int(int(exp(-x^2-z^2),'x',z-zlim,zlim-z),'z',-inf,zlim)
 
ans =
 
int(-pi^(1/2)*erf(z - zlim)*exp(-z^2), z == -Inf..zlim)

Not a solution that doesn't require an integration - but at least something like a solution...

Answer 1

[Bjorn Gustavsson]

Your problem might be that the symbolic toolbox have no way of knowing that your correlation is between -1 and 1. If you try the trivial calculation for the case where there is no correlation it is unproblematically fast. One way forward might be to transform your coordinates to one where there is no correlation. This would have to be combined with a change in the integration-limits, but by a simple sketch this seems "reasonably" manageable to do - the rectangular region in your original variables ought to become a triangular region in the transformed variables.

HTH

[Donghwan Hyun]

Thanks a lot