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.
回答（2 个）
Matt Tearle
2011-3-23
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 :)
Walter Roberson
2011-3-23
0 个投票
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.
4 个评论
Walter Roberson
2011-3-23
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
2011-3-23
I got a little less than 1/16 numerically -- about 0.055.
Walter Roberson
2011-3-23
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
2011-3-23
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)
类别
在 帮助中心 和 File Exchange 中查找有关 Hypergeometric Distribution 的更多信息
2011-3-23
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
