[数学]二维对数正态分布的概率分布，期望，方差和相关系数

最近遇到了一个联合对数正态分布的相关系数的问题，搜遍全网无果，索性自己动手。本文借鉴了这个知乎回答

首先我们有二维正态分布：
$$X,Y\sim \mathbf{BVN}({\mu }_x,{\mu }_y,{\sigma }_x^2,{\sigma }_y^2,{\rho }_{xy})$$

取对数之后我们会得到二维对数正态分布的概率密度函数。只写了第一象限的函数表达式，其他地方都是0。
$$f(x,y)=\frac{1}{2\pi \sqrt{1-{\rho }_{xy}^2}{\sigma }_x{\sigma }_{y}xy}\exp \left[\frac{-1}{2(1-{\rho }_{xy}^2)}\left(\frac{(\ln x-{\mu }_x{)}^2}{{\sigma }_x^2}-\frac{2{\rho }_{xy}(\ln x-{\mu }_x)(\ln y-{\mu }_y)}{{\sigma }_x{\sigma }_y}+\frac{(\ln y-{\mu }_y{)}^2}{{\sigma }_y^2}\right)\right]$$

引用链接里有边缘分布（一维情况下）的期望和方差的推导过程，这里只写结论：
$$\begin{gathered} E(X)=\exp ({\mu }_x+\frac{{\sigma }_x^2}{2}) \\ D(X)=\exp (2{\mu }_x+{\sigma }_x^2)(\exp ({\sigma }_x^2)-1) \end{gathered}$$

接下来想算相关系数。首先我们有相关系数的公式：
$$\rho =\frac{COV(X,Y)}{\sqrt{D(X)D(Y)}}=\frac{E(XY)-E(X)E(Y)}{\sqrt{D(X)D(Y)}}$$

关键一步是计算$E(XY)$
$$E(XY)={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }xyf(x,y)\mathbf{d}x\mathbf{d}y$$

代入$f(x,y)$
$$E(XY)={\int }_0^{+\infty }{\int }_0^{+\infty }\frac{1}{2\pi \sqrt{1-{\rho }_{xy}^2}{\sigma }_x{\sigma }_y}\exp \left[\frac{-1}{2(1-{\rho }_{xy}^2)}\left(\frac{(\ln x-{\mu }_x{)}^2}{{\sigma }_x^2}-\frac{2{\rho }_{xy}(\ln x-{\mu }_x)(\ln y-{\mu }_y)}{{\sigma }_x{\sigma }_y}+\frac{(\ln y-{\mu }_y{)}^2}{{\sigma }_y^2}\right)\right]\mathbf{d}x\mathbf{d}y$$

作变换（"简单的"二次型标准化）
$$u=\frac{\ln x-{\mu }_x}{{\sigma }_x}-({\rho }_{xy}{\sigma }_y+{\sigma }_x),v=\frac{\ln y-{\mu }_y}{{\sigma }_y}-({\rho }_{xy}{\sigma }_x+{\sigma }_y)$$

逆变换及其微分
$$\begin{gathered} x=\exp ({\sigma }_{x}u+{\rho }_{xy}{\sigma }_x{\sigma }_y+{\sigma }_x^2+{\mu }_x), \\ y=\exp ({\sigma }_{y}u+{\rho }_{xy}{\sigma }_x{\sigma }_y+{\sigma }_y^2+{\mu }_y), \\ \mathbf{d}x={\sigma }_x\exp ({\sigma }_{x}u+{\rho }_{xy}{\sigma }_x{\sigma }_y+{\sigma }_x^2+{\mu }_x)\mathbf{d}u, \\ \mathbf{d}y={\sigma }_y\exp ({\sigma }_{y}u+{\rho }_{xy}{\sigma }_x{\sigma }_y+{\sigma }_y^2+{\mu }_y)\mathbf{d}u. \end{gathered}$$

代入$E(XY)$得（节省空间不写积分上下限了，都是无穷）
$$\begin{gathered} E(XY)=\frac{1}{2\pi \sqrt{1-{\rho }_{xy}^2}}\iint \exp \left[\frac{-1}{2(1-{\rho }_{xy}^2)}(u+{\rho }_{xy}{\sigma }_y+{\sigma }_x{)}^2-2{\rho }_{xy}(u+{\rho }_{xy}{\sigma }_y+{\sigma }_x)(v+{\rho }_{xy}{\sigma }_x+{\sigma }_y) \\ +(v+{\rho }_{xy}{\sigma }_x+{\sigma }_y{)}^2+{\sigma }_{x}u+{\sigma }_{y}v+2{\rho }_{xy}{\sigma }_x{\sigma }_y+{\sigma }_x^2+{\sigma }_y^2+{\mu }_x+{\mu }_y\right]\mathbf{d}u\mathbf{d}v \end{gathered}$$

化简得到
$$E(XY)=\exp ({\mu }_x+{\mu }_y+\frac{1}{2}({\sigma }_x^2+2{\rho }_{xy}{\sigma }_x{\sigma }_y+{\sigma }_y^2))\frac{1}{2\pi (1-{\rho }_{xy}^2)}\iint \exp \left[\frac{-1}{2(1-{\rho }_{xy}^2)}(u^2-2{\rho }_{xy}uv+v^2)\right]\mathbf{d}u\mathbf{d}v$$

指数项右边是一个正态分布概率密度的积分，因此等于1，于是得到了一个很简单的形式
$$E(XY)=\exp ({\mu }_x+{\mu }_y+\frac{1}{2}({\sigma }_x^2+2{\rho }_{xy}{\sigma }_x{\sigma }_y+{\sigma }_y^2))$$

然后我们把$E(XY)$，$E(X)$，$E(Y)$，$D(X)$，$D(Y)$代入相关系数公式化简得

$$\rho =\frac{\exp \left({\rho }_{xy}{\sigma }_x{\sigma }_y\right)-1}{\sqrt{(\exp ({\sigma }_x^2)-1)(\exp ({\sigma }_y^2)-1)}}$$

但是这个相关系数的结果有个很奇怪的性质，困扰了我一天，那就是当${\sigma }_x\ne {\sigma }_y$的时候$\rho$取不到$[-1,1]$，我用数字帝国画了个${\sigma }_x=1，{\sigma }_y=2$时的草图，长这样：
然后就怀疑我哪里做错了，后来想着还是拿matlab数值计算一下。代码如下：

rho = 0.99;
sigma_x = 2;
sigma_y = 1;
mu_x = 1;
mu_y = 1;
%ff = @(x,y)(exp(-((((log(x)-mu_x).^2./sigma_x.^2)-(2.*rho.*(log(x)-mu_x).*(log(y)-mu_y)./(sigma_x.*sigma_y))+((log(y)-mu_y).^2)./sigma_y.^2)./(2.*(1-rho.^2))))./(2*sigma_x*sigma_y.*pi.*sqrt(1-rho.^2).*x.*y));原始函数
fexy = @(x, y)(exp(-((((log(x)-mu_x).^2./sigma_x.^2)-(2.*rho.*(log(x)-mu_x).*(log(y)-mu_y)./(sigma_x.*sigma_y))+((log(y)-mu_y).^2)./sigma_y.^2)./(2.*(1-rho.^2))))./(2*sigma_x*sigma_y.*pi.*sqrt(1-rho.^2)));
exy = integral2(fexy,0,inf,0,inf,'Method','iterated','AbsTol',0,'RelTol',1e-10);
exey = exp(mu_x+mu_y+sigma_x^2/2+sigma_y^2/2);
corr = (exy-exey)/(exey*sqrt((exp(sigma_x^2)-1)*(exp(sigma_y^2)-1)));

结果是0.6505，和图像相符，也就是说二维对数正态分布的相关系数取值范围确实不总是$[-1,1]$。
再附一个画二维正态和二维对数正态概率分布的代码：

X1=[0.01:0.01:3];
Y1=[0.01:0.01:3];
[x,y]=meshgrid(X1,Y1);
rho = 0.5;
sigma_x = 1;
sigma_y = 1;
mu_x = 1;
mu_y = 1;
BVLN=(exp(-((((log(x)-mu_x).^2./sigma_x.^2)-(2.*rho.*(log(x)-mu_x).*(log(y)-mu_y)./(sigma_x.*sigma_y))+((log(y)-mu_y).^2)./sigma_y.^2)./(2.*(1-rho.^2))))./(2*sigma_x*sigma_y.*pi.*sqrt(1-rho.^2).*x.*y));
BVN=(exp(-((((x-mu_x).^2./sigma_x.^2)-(2.*rho.*(x-mu_x).*(y-mu_y)./(sigma_x.*sigma_y))+((y-mu_y).^2)./sigma_y.^2)./(2.*(1-rho.^2))))./(2*sigma_x*sigma_y.*pi.*sqrt(1-rho.^2)));
subplot(1,2,1);surf(x,y,BVLN);
subplot(1,2,2);surf(x,y,BVN);

画出来是这种感觉：

断断续续算了四天（主要是开始时不知道如何做变换），算的心态爆炸，给个免费的赞再走吧。