solving non-linear function with double integral

Would any one please help me on how to code this up?

I want to solve 'R' for

int(int(1/(2*pi*sqrt(1-R^2))*exp(-(1/(2*(1-R^2))*(x^2-2*R*x*y+y^2))),x, -inf, -2.7450), y, -inf, -2.7450) - 0.0000193 = 0

basically, it is a bivariate normal distribution with mean at [0, 0] but the correlation factor R is the unknown. I tried to use mvncdf function but it seems to me that it doesn't take unknown correlation factor.

Thanks.

Answer 1

[Matt Tearle]

How's this:

%%Find R
g1 = @(x,y,R) 1./(2*pi*sqrt(1-R.^2)).*exp(-(1./(2*(1-R.^2)).*(x.^2-2*R.*x.*y+y.^2)));
g2 = @(R) dblquad(@(x,y) g1(x,y,R),-100,-2.7,-100,-2.7,1e-8)- 0.0000193;
g3 = @(R) quad2d(@(x,y) g1(x,y,R),-100,-2.7,-100,-2.7,'AbsTol',1e-8,'RelTol',1e-5)- 0.0000193;

R = fzero(g2,0)
R = fzero(g3,0)

%%Sanity check
R = linspace(-0.1,0.2);
I1 = zeros(size(R));
I2 = zeros(size(R));
for k=1:length(I1)
I1(k) = g2(R(k));
I2(k) = g3(R(k));
end
plot(R,I1,R,I2)
grid

Obviously, this is more than you asked for. Just trying to show that quad2d and dblquad give the same answers, given sufficiently tight error tolerances. All you really need is lines 1, 3, and 5.

Note the use of -100 in place of -Inf. Replace with -100000 if you like :)

Answer 2

[Walter Roberson]

Symbolically, the double integral converts to a single integral over y, of an expression which is the limit as x approaches negative infinity. x appears only once in the expression, in the numerator of a fraction in an erf() call. The fraction is singular at R=-1 and R=1, but continuous in (-1,1) and the denominator has a consistent sign over that range. You end up doing a double or triple negation of the infinity in the numerator. Assuming R in (-1,1) the limit of the overall expression then becomes the overall expression with erfc(-infinity) or erfc(infinity) in place of that term [not sure which, do the sign analysis!] which gives you a -1 or +1 for that term. After that substitution the expression converts in to a single integral over y.

The single integral can be simplified to reduce the complexity of th erfc() call that remains, especially if you convert to hypergeometric.

Unfortunately it looks likely the single integral would have to be evaluated numerically, and doing that for even a single R value appears to be slow. You could use fzero() on it in theory, but Matt's approach might be much much faster in practice.

[Walter Roberson]

The extreme slowness appears to be below R=-0.7; at that value, the integral itself is very close to 0, below that it appears to be non-integrable. From -0.7 upwards, individual evaluations are not too bad. The overall solution is somewhere between 1/16 and 1/8

[Matt Tearle]

I got a little less than 1/16 numerically -- about 0.055.

[Walter Roberson]

If you substitute in 1/16 before the integration then the overall expression fairly quickly evaluates to -0.342670685e-5 but at 1/8 evaluates to a positive number, so the root is between the two.

[Walter Roberson]

0.08633825363 I make it. The Maple expression used was

fsolve(int(int(exp(-(x^2-2*R*x*y+y^2)/(2*(1-R^2)))/(2*Pi*sqrt(1-R^2)), x = -infinity .. -2.7450), y = -infinity .. -2.7450)-0.193e-4, R = 1/16 .. 1/8)