matlab - 双变量对数正态密度函数的积分

Question

rho = 0.8;

ff = @(x, y) (exp(-(((log(x)-10).^2 - 2.* rho .* (log(x)-10) .* (log(y)-10)+(log(y)-10).^2)./(2 .* (1-rho.^2))))./(2.*pi.*sqrt(1-rho.^2).*x.*y));

syms x y
vpaintegral(vpaintegral(ff, x, [0 inf]), y, [0 inf])

为什么上面的二元对数正态密度函数在Matlab中的积分不是1？

注意:此对数正态二元随机变量的对数变换是一个均值 (10, 10) 和协方差矩阵 (1, rho, rho, 1) 的二元正态随机变量。

Answer 1

使用 integral2 :

我们得到 0.9994，

% MATLAB R2019a
rho = 0.8;
ff = @(x, y) (exp(-(((log(x)-10).^2-2.*rho.*(log(x)-10).*(log(y)-10)+(log(y)-10).^2)./(2.*(1-rho.^2))))./(2.*pi.*sqrt(1-rho.^2).*x.*y));

area = integral2(ff,0,inf,0,inf)     % area = 0.9994

但调整容差可以得到所需的结果。

area = integral2(ff,0,inf,0,inf,'Method','iterated','AbsTol',0,'RelTol',1e-10)

ans = 1.0000

format long
area

ans = 0.999999999999998

不算太破烂。

---

使用 vpaintegral 来自符号工具箱:

您还可以adjust the tolerance对于 vpaintegral .

使用 Relative Error Tolerance 1e-4 完成了工作。该参数极大地影响计算时间。

syms x y
area = vpaintegral(vpaintegral(ff, x, [0 inf],'RelTol', 1e-4, 'AbsTol', 0), y, [0 inf],'RelTol', 1e-4, 'AbsTol', 0)

area = 1.0