GitModel-动手学概率论_02（python）

1 动手学概率_2

github 上pdf版本及ipynb版本：https://github.com/cx-333/Math-Modeling

1.5 多维随机变量及其联合分布、边际分布、条件分布

• 多维随机变量：若随机变量 $X_1(\omega ),X_2(\omega ),\cdots ,X_n(\omega )$ 定义在同一个基本空间 Ω = { ω } \Omega=\{\omega\} Ω={ω} 上， 则称
  $$\mathbf{X}(\omega )=\left(X_1(\omega ),X_2(\omega ),\cdots ,X_n(\omega )\right)$$
  是一个多维随机变量，也称为n维随机向量。

• 多维随机变量联合分布函数：设 $X=\left(X_1,X_2,\cdots ,X_n\right)$ 是 $n$ 维随机变量， 对任意 $n$ 个实数 $x_1,x_2,\cdots ,x_n$ 所组成的 $n$ 个事件 $X_{1} \leqslant x_{1},X_{2} \leqslant x_{2} , \cdots, X_{n} \leqslant x_{n} $ 同时发生的概率
  $$F\left(x_1,x_2,\cdots ,x_n\right)=P\left(X_1⩽x_1,X_2⩽x_2,\cdots ,X_n⩽x_n\right)$$
  称为 $n$ 维随机变量 $\mathbf{X}$ 的联合分布函数。

• 多维随机变量联合概率密度函数：设$n$维随机变量 $X=\left(X_1,X_2,\cdots ,X_n\right)$ 的分布函数为 $F(x_1,x_2,\cdots ,x_n)$ 。假如各分量 $x_1,x_2,\cdots ,x_n$ 都是一维连续随机变量，并存在定义在空间上的非负函数 $p(x_1,x_2,\cdots ,x_n)$，使得
  $$F(x,y)={\int }_{-\infty }^{x_1}{\int }_{-\infty }^{x_2}\cdots {\int }_{-\infty }^{x_n}p(x_1,x_2,\cdots ,x_n)d{x}_{1}d{x}_2\cdots d{x}_n$$
  则称 $X=\left(X_1,X_2,\cdots ,X_n\right)$ 为$n$维连续随机变量，$p(x_1,x_2,\cdots ,x_n)$ 称为 $X=\left(X_1,X_2,\cdots ,X_n\right)$ 的联合概率密度函数， 或简称联合密度。

• 特例：二维随机变量
  二维随机变量的分布函数：设$(X,Y)$是二维随机变量，对于任意实数$x,y$，二元函数：
  F ( x , y ) = P { X ≤ x , Y ≤ y } F(x, y)=P\{X\le x, Y \le y\} F(x,y)=P{X≤x,Y≤y}
  称为二维随机变量$(X,Y)$的分布函数，或称为随机变量$(X,Y)$的联合分布函数。

分布函数的性质：

• 1. $F(x,y)$是$x,y$的不减函数
• 2. $0\le F(x,y)\le 1$，且$$\begin{array}{rl}& 对于任意固定的y,F(-\infty ,y)=0\\ & 对于任意固定的x,F(x,-\infty )=0\\ & F(-\infty ,-\infty )=0,F(\infty ,\infty )=1\end{array}$$
• 3. $F(x+0,y)=F(x,y),F(x,y+0)=F(x,y)$，即$F(x,y)$关于$x,y$分别右连续
• 4. 对于任意的$(x_1,y_1),(x_2,y_2),x_1<x_2,y_1<y_2$，有不等式$$F(x_2,y_2)-F(x_2,y_1)+F(x_1,y_1)-F(x_1,y_2)\ge 0$$恒成立

  二维随机变量的概率密度函数：对于二维随机变量$(X,Y)$的分布函数$F(x,y)$，如果存在非负可积函数$f(x,y)$，使对于任意$x,y$有
$$F(x,y)={\int }_{-\infty }^y{\int }_{-\infty }^{x}f(u,v)dudv$$
则称$(X,Y)$是连续型的二维随机变量，函数$f(x,y)$称为二维随机变量$(X,Y)$的概率密度，或称为随机变量$X和Y$的联合概率密度。

联合概率密度$f(x,y)$的性质：

• 1. $f(x,y)ge0$
• 2. ${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)dxdy=1$
• 3. 若$f(x,y)$在点$(x,y)$连续，则有$$\frac{{\partial }^{2}F(x,y)}{\partial x\partial y}=f(x,y)$$
• 4. 设$G$是$xOy$平面上的区域，点$(X,Y)$落在$G$内的概率为： P { ( X , Y ) ∈ G } = ∬ G f ( x , y ) d x d y P\{(X, Y) \in G\} = \iint\limits_{G}f(x, y)dxdy P{(X,Y)∈G}=G∬​f(x,y)dxdy

python代码（绘制二维概率密度函数图像）代码参考

import numpy as np
from scipy import stats
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt

x, y = np.mgrid[-5:5:0.01, -5:5:0.01]
pos = np.dstack((x, y))
rv = stats.multivariate_normal([0.5, -0.2], [[2.0, 0.3], [0.3, 0.5]])
z = rv.pdf(pos)
# 曲面图
plt.figure('Surface', facecolor='lightgray', figsize=(12, 8))
ax = plt.axes(projection='3d')
ax.set_xlabel('x', fontsize=14)
ax.set_ylabel('y', fontsize=14)
ax.set_zlabel('f(x, y)', fontsize = 14)
ax.plot_surface(x, y, z, rstride=50, cstride=50, cmap='jet')
plt.show()
# 等高线图
fig2 = plt.figure(figsize=(8,6))
ax2 = fig2.add_subplot(111)
ax2.contourf(x, y, rv.pdf(pos), rstride=50, cstride=50, cmap='jet')
plt.show()

C:\Users\13541\AppData\Local\Temp\ipykernel_11684\2054532681.py:21: UserWarning: The following kwargs were not used by contour: 'rstride', 'cstride'
  ax2.contourf(x, y, rv.pdf(pos), rstride=50, cstride=50, cmap='jet')

• 边际（边缘）分布：
  二维随机变量$(X,Y)$作为一个整体，具有分布函数$F(x,y)$，$X和Y$都是随机变量，他们也有自己的分布函数，分别记为$F_X(x),F_Y(y)$，依次为随机变量$(X,Y)$关于$X$和关于$Y$的边缘分布函数 F X ( x ) = P { X ≤ x } = P { X ≤ x , Y < ∞ } = F ( x , ∞ ) F_{X}(x)=P\{X\le x\}=P\{X\le x, Y < \infty\}=F(x, \infty) FX​(x)=P{X≤x}=P{X≤x,Y<∞}=F(x,∞)即$F_X(x)=F(x,\infty )$，同理$$F_Y(y)=F(\infty ,y)$$
  对于连续型随机变量$(X,Y)$，设它的概率密度函数为$f(x,y)$，由于 $$F_X(x)=F(x,\infty )={\int }_{-\infty }^x[{\int }_{-\infty }^{\infty }f(x,y)dy]dx$$则随机变量$X$的概率密度函数为：$f_X(x)={\int }_{-\infty }^{\infty }f(x,y)dy$
  同理，随机变量$Y$的概率密度函数为：$f_Y(y)={\int }_{-\infty }^{\infty }f(x,y)dx$

python代码（已知联合分布求边缘分布）

from sympy import *

x, y = symbols('x, y')
fxy = 1
pxy = Piecewise((fxy, And(x > 0, x < 1, y > -x, y < x)), (0, True))
fXx = integrate(pxy, (y, -oo, oo))
pprint("关于随机变量X的边缘密度函数为：{}".format(fXx))
fYy = integrate(pxy, (x, -oo, oo))
pprint("关于随机变量Y的边缘密度函数为：{}".format(fYy))

关于随机变量X的边缘密度函数为：Piecewise((x + Max(-x, x), (x > 0) & (x < 1)), (0, True))
关于随机变量Y的边缘密度函数为：-Max(0, -y, y) + Max(1, -y, y)

  边际（边缘）分布列：在二维离散随机变量 $(X,Y)$ 的联合分布列 { P ( X = x i , Y = y j ) } \left\{P\left(X=x_{i}, Y=y_{j}\right)\right\} {P(X=xi​,Y=yj​)} 中， 对 $j$ 求和所得的分布列$$\sum _{j=1}^{\infty }P\left(X=x_i,Y=y_j\right)=P\left(X=x_i\right),i=1,2,\cdots$$
  被称为 $X$ 的边际分布列。 类似地， 对 $i$ 求和所得的分布列 $$\sum _{i=1}^{\infty }P\left(X=x_i,Y=y_j\right)=P\left(Y=y_j\right),j=1,2,\cdots$$
  被称为 $Y$ 的边际分布列。

• 条件分布：由条件概率引出条件概率分布。
  设$(X,Y)$是二维离散型随机变量，其分布律为 P { X = x i , Y = y j } = p i j , i , j = 1 , 2 , ⋯   . P\{X=x_{i}, Y=y_{j}\}=p_{ij}, i,j=1,2,\cdots . P{X=xi​,Y=yj​}=pij​,i,j=1,2,⋯. $(X,Y)$关于$X$和关于$Y$的边缘分布分别为
  P { X = x i } = p i = ∑ j = 1 ∞ p i j , i = 1 , 2 , ⋯   . P { Y = y j } = p j = ∑ i = 1 ∞ p i j , j = 1 , 2 , ⋯   . \begin{aligned}& P\{X=x_{i}\}=p_{i}=\sum_{j=1}^{\infty}p_{ij}, i=1,2,\cdots .\\&P\{Y=y_{j}\}=p_{j}=\sum_{i=1}^{\infty}p_{ij}, j=1,2,\cdots . \end{aligned} ​P{X=xi​}=pi​=j=1∑∞​pij​,i=1,2,⋯.P{Y=yj​}=pj​=i=1∑∞​pij​,j=1,2,⋯.​
  设$P_{\cdot j}>0$，考虑在事件 { Y = y j } \{Y=y_{j}\} {Y=yj​}已发生的条件下事件 { X = x i } \{X=x_{i}\} {X=xi​}发生的概率，由条件概率公式得 P { X = x i ∣ Y = y j } = P { X = x i , Y = y j } P { Y = y j } = p i j p ⋅ j , i = 1 , 2 , ⋯   . P\{X=x_{i}|Y=y_{j}\} = \frac{P\{X=x_{i}, Y=y_{j}\}}{P\{Y=y_{j}\}} = \frac{p_{ij}}{p_{\cdot j}}, i=1,2,\cdots . P{X=xi​∣Y=yj​}=P{Y=yj​}P{X=xi​,Y=yj​}​=p⋅j​pij​​,i=1,2,⋯.上式称为在$Y=y_j$条件下随机变量$X$得条件分布率，同理 P { Y = y j ∣ X = x i } = P { Y = y j , X = x i } P { X = x i } = p i j p i ⋅ , i = 1 , 2 , ⋯   . P\{Y=y_{j}|X=x_{i}\} = \frac{P\{Y=y_{j}, X=x_{i}\}}{P\{X=x_{i}\}} = \frac{p_{ij}}{p_{i \cdot}}, i=1,2,\cdots . P{Y=yj​∣X=xi​}=P{X=xi​}P{Y=yj​,X=xi​}​=pi⋅​pij​​,i=1,2,⋯.称为在$X=x_i$条件下随机变量$Y$得条件分布率

条件分布得性质：

• 1. P { X = x i ∣ Y = y j } ≥ 0 P\{X=x_{i}|Y=y_{j}\} \ge 0 P{X=xi​∣Y=yj​}≥0
• 2. ∑ i = 1 ∞ P { X = x i ∣ Y = y j } = ∑ i = 1 ∞ p i j p ⋅ j = 1 p ⋅ j ∑ i = 1 ∞ p i j = p ⋅ j p ⋅ j = 1 \sum_{i=1}^{\infty}P\{X=x_{i}|Y=y_{j}\} = \sum_{i=1}^{\infty}\frac{p_{ij}}{p_{\cdot j}} =\frac{1} {p_{\cdot j}}\sum_{i=1}^{\infty}p_{ij} = \frac{p_{\cdot j}}{p_{\cdot j}} = 1 i=1∑∞​P{X=xi​∣Y=yj​}=i=1∑∞​p⋅j​pij​​=p⋅j​1​i=1∑∞​pij​=p⋅j​p⋅j​​=1

  连续型随机变量条件分布：
设二维随机变量$(X,Y)$得概率密度为$f(x,y)$，$(X,Y)$关于$Y$得边缘概率密度为$f_Y(y)$. 若对于固定得$y,f_Y(y)>0$，则称$\frac{f(x,y)}{f_Y(y)}$为在$Y=y$条件下$X$的条件概率密度函数，记为
$$f_{X|Y}(x|y)=\frac{f(x,y)}{f_Y(y)}$$
称${\int }_{-\infty }^x{f}_{X|Y}(x|y)dx={\int }_{-\infty }^x\frac{f(x,y)}{f_Y(y)}dx$为在$Y=y$条件下$X$的条件分布函数，记为 P { X ≤ x ∣ Y = y } 或 F X ∣ Y ( x ∣ y ) P\{X\le x|Y =y\}或F_{X|Y}(x|y) P{X≤x∣Y=y}或FX∣Y​(x∣y)，即
F X ∣ Y ( x ∣ y ) = P { X ≤ x ∣ Y = y } = ∫ − ∞ x f ( x , y ) f Y ( y ) d x . F_{X|Y}(x|y)=P\{X\le x|Y =y\}=\int_{-\infty}^{x}\frac{f(x, y)}{f_{Y}(y)}dx. FX∣Y​(x∣y)=P{X≤x∣Y=y}=∫−∞x​fY​(y)f(x,y)​dx.

连续随机变量的贝叶斯公式和全概率公式：

• 1. 全概率公式$$\begin{array}{rl}& p_Y(y)={\int }_{-\infty }^{\infty }{p}_X(x)p(y\mid x)\mathrm{d}x,\\ & p_{\chi }(x)={\int }_{-\infty }^{\infty }{p}_Y(y)p(x\mid y)\mathrm{d}y.\end{array}$$
• 2. 贝叶斯公式$$\begin{array}{rl}& p(x\mid y)=\frac{p_X(x)p(y\mid x)}{{\int }_{-\infty }^{\infty }{p}_X(x)p(y\mid x)\mathrm{d}x},\\ & p(y\mid x)=\frac{p_Y(y)p(x\mid y)}{{\int }_{-\infty }^{\infty }{p}_Y(y)p(x\mid y)\mathrm{d}y}.\end{array}$$

python代码(求解边际分布列)：设在一段时间内进人某一商店的顾客人数 $X$ 服从泊松分布 $P(\lambda )$， 每个顾客购买某种物品的概率为 $p$, 并且各个顾客是否购买该种物品相互独立， 求进入商店的顾客购买这种物品的人数 $Y$ 的分布列。

from sympy import *
from sympy.abc import lamda, k, m, x, y, p    # 代替symbols("lamda, k, m, x, y") 
Pxk = (lamda ** k * exp(- lamda)) / factorial(k)
print("进入商店人数(k)的概率分布P(X=k) = {}".format(Pxk))
Pyx = (factorial(k) / (factorial(m) * factorial(k-m))) * (p ** m) * ((1 - p) ** (k - m))
print("在进入商店人数（k）确定的条件下，客户买某种商品人数的条件分布P(Y=m|X =k)={}".format(Pyx))
f = Pxk * Pyx
Py = summation(f, (k, m, oo))
print("进入商店的顾客购买这种商品人数Y的分布列P(Y=m)={}".format(Py))
Py

进入商店人数(k)的概率分布P(X=k) = lamda**k*exp(-lamda)/factorial(k)
在进入商店人数（k）确定的条件下，客户买某种商品人数的条件分布P(Y=m|X =k)=p**m*(1 - p)**(k - m)*factorial(k)/(factorial(m)*factorial(k - m))
进入商店的顾客购买这种商品人数Y的分布列P(Y=m)=lamda**m*p**m*exp(-lamda)*exp(-lamda*(p - 1))/factorial(m)

$\frac{{\lambda }^m{p}^m{e}^{-\lambda }{e}^{-\lambda \left(p-1\right)}}{m!}$

1.6 多维随机变量的数字特征：期望向量、协方差与协方差矩阵、相关系数与相关系数矩阵

• 期望向量
  记 $n$ 维随机向量为 $\mathbf{X}={\left(X_1,X_2,\cdots ,X_n\right)}^T$, 若其每个分量的数学期望都存在， 则称
  $$E(\mathbf{X})={\left(E\left(X_1\right),E\left(X_2\right),\cdots ,E\left(X_n\right)\right)}^T$$
  为 $n$ 维随机向量 $\mathbf{X}$ 的数学期望向量（一般为列向量）， 简称为 $\mathbf{X}$ 的数学期望。

• 协方差与协方差矩阵
    协方差：$\mathrm{Cov}(X,Y)=E[(X-E(X))(Y-E(Y))]$ ，衡量的是两个随机变量之间的相互关联的程度

1. 当 $\mathrm{Cov}(X,Y)>0$ 时， 称 $X$ 与 $Y$ 正相关， 这时两个偏差 $(X-E(X))$ 与 $(Y-E(Y))$ 有同时增加或同时减少的倾向。 由于 $E(X)$ 与 $E(Y)$ 都是常数， 故等价于 $X$ 与 $Y$ 有同时增加或同时减少的倾向。
2. 当 $\mathrm{Cov}(X,Y)<0$ 时， 称 $X$ 与 $Y$ 负相关, 这时有 $X$ 增加而 $Y$ 减少的倾向， 或有 $Y$ 增加而 $X$ 减少的倾向。
3. 当 $\mathrm{Cov}(X,Y)=0$ 时，称 $X$ 与 $Y$ 不相关。 这时可能由两类情况导致:一类是 $X$ 与 $Y$ 的取值毫无关联， 另一类是 $X$ 与 $Y$ 间存有某种非线性关系。

协方差$\mathrm{Cov}(X,Y)$的性质：

• 1. $\mathrm{Cov}(X,Y)=E(XY)-E(X)E(Y)$
• 2. 若随机变量 $X$ 与 $Y$ 相互独立， 则 $\mathrm{Cov}(X,Y)=0$， 反之不成立。
• 3. （最重要）对任意二维随机变量 $(X,Y)$， 有$$\mathrm{Var}(X\pm Y)=\mathrm{Var}(X)+\mathrm{Var}(Y)\pm 2\mathrm{Cov}(X,Y)$$ 该性质表明: 在 $X$ 与 $Y$ 相关的场合,和的方差不等于方差的和。 $X$ 与 $Y$ 的正相关会增加和的方差,负相关会减少和的方差，而在 $X$ 与 $Y$ 不相关的场合，和的方差等于方差的和，即：若 $X$ 与 $Y$ 不相关， 则 $\mathrm{Var}(X\pm Y)=\mathrm{Var}(X)+\mathrm{Var}(Y)$。
• 4. 协方差 $\mathrm{Cov}(X,Y)$ 的计算与 $X,Y$ 的次序无关， 即 $$\mathrm{Cov}(X,Y)=\mathrm{Cov}(Y,X)$$
• 5. 任意随机变量 $X$ 与常数 $a$ 的协方差为零，即 $$\mathrm{Cov}(X,a)=0$$
• 6. 对任意常数 $a,b$， 有 $$\mathrm{Cov}(aX,bY)=ab\mathrm{Cov}(X,Y)$$
• 7. 设 $X,Y,Z$ 是任意三个随机变量,则 $$\mathrm{Cov}(X+Y,Z)=\mathrm{Cov}(X,Z)+\mathrm{Cov}(Y,Z)$$

🔥例题：设随机变量$(X,Y)$既有概率密度
$$f(x,y)=\{\begin{array}{rlr}& \frac{1}{8}(x+y),& 0\le x\le 2,0\le y\le 2\\ & 0,& 其它\end{array}$$
求$E(X),E(Y),E(XY),Cov(X,Y),D(X+Y)$

✨解：
$$\begin{array}{rl}& E(X)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }xf(x,y)dxdy\\ & E(Y)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }yf(x,y)dxdy\\ & E(XY)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }xyf(x,y)dxdy\\ & Cov(X,Y)=E(XY)-E(X)E(Y)\\ & D(X+Y)=D(X)+D(Y)+2Cov(X,Y)=E(X^2)-[E(X){]}^2+E(Y^2)-[E(Y){]}^2+2Cov(X,Y)\end{array}$$

python代码(求解例题)

from sympy import *
from sympy.abc import x, y 
f = (1 / 8) * (x + y)
fxy = Piecewise((f, And(x >= 0, x <= 2, y >= 0, y <= 2)), (0, True))
Ex = integrate(x * fxy, (x, -oo, oo), (y, -oo, oo))
print("E(X) = {}".format(Ex))
Ey = integrate(y * fxy, (x, -oo, oo), (y, -oo, oo))
print("E(Y) = {}".format(Ey))
Exy = integrate(x * y * fxy, (x, -oo, oo), (y, -oo, oo))
print("E(XY) = {}".format(Exy))
cov_xy = Exy - Ex * Ey 
print("Cov(X, Y) = {}".format(cov_xy))
Ex_2 = integrate(x ** 2 * fxy, (x, -oo, oo), (y, -oo, oo))
Ey_2 = integrate(y ** 2 * fxy, (x, -oo, oo), (y, -oo, oo))
Dx_y = Ex_2 - (Ex ** 2) + Ey_2 - (Ey ** 2) + 2 * cov_xy
print("D(X + Y) = {}".format(Dx_y))

E(X) = 1.16666666666667
E(Y) = 1.16666666666667
E(XY) = 1.33333333333333
Cov(X, Y) = -0.0277777777777775
D(X + Y) = 0.555555555555557

  协方差矩阵：设$n$ 维随机向量为 $\mathbf{X}={\left(X_1,X_2,\cdots ,X_n\right)}^{\prime }$的期望向量为：
$$E(\mathbf{X})={\left(E\left(X_1\right),E\left(X_2\right),\cdots ,E\left(X_n\right)\right)}^T$$
则
$$\begin{array}{rl}& E\left[(\mathbf{X}-E(\mathbf{X}))(\mathbf{X}-E(\mathbf{X}){)}^T\right]\\ =& \left(\begin{array}{cccc}\mathrm{Var}\left(X_1\right)& \mathrm{Cov}\left(X_1,X_2\right)& \cdots & \mathrm{Cov}\left(X_1,X_n\right)\\ \mathrm{Cov}\left(X_2,X_1\right)& \mathrm{Var}\left(X_2\right)& \cdots & \mathrm{Cov}\left(X_2,X_n\right)\\ ⋮& ⋮& & ⋮\\ \mathrm{Cov}\left(X_n,X_1\right)& \mathrm{Cov}\left(X_n,X_2\right)& \cdots & \mathrm{Var}\left(X_n\right)\end{array}\right)\end{array}$$
称为该随机向量的方差-协方差矩阵，简称协方差阵，记为 $\mathrm{Cov}(\mathbf{X})$。

🦊注：$n$ 维随机向量的协方差矩阵 $\mathrm{Cov}(\mathbf{X})={\left(\mathrm{Cov}\left(X_i,X_j\right)\right)}_{n\times n}$ 是一个对称的非负定矩阵。

python代码（求上一例题的协方差矩阵）

from sympy import *
from sympy.abc import x, y 
f = (1 / 8) * (x + y)
fxy = Piecewise((f, And(x >= 0, x <= 2, y >= 0, y <= 2)), (0, True))
Ex = integrate(x * fxy, (x, -oo, oo), (y, -oo, oo))
Ey = integrate(y * fxy, (x, -oo, oo), (y, -oo, oo))
Exy = integrate(x * y * fxy, (x, -oo, oo), (y, -oo, oo))

cov_xy = Exy - Ex * Ey 
var_x = Ex_2 - (Ex ** 2)
var_y = Ey_2 - (Ey ** 2)
Matrix([[var_x, cov_xy], [cov_xy, var_y]])

$\left[\begin{array}{cc}0.305555555555556& -0.0277777777777775\\ -0.0277777777777775& 0.305555555555556\end{array}\right]$

• 相关系数与相关系数矩阵
    相关系数：设 $(X,Y)$ 是一个二维随机变量， 且 $\mathrm{Var}(X)={\sigma }_X^2>0,\mathrm{Var}(Y)={\sigma }_Y^2>0$.则称
  $$\mathrm{Corr}(X,Y)=\frac{\mathrm{Cov}(X,Y)}{\sqrt{\mathrm{Var}(X)}\sqrt{\mathrm{Var}(Y)}}=\frac{\mathrm{Cov}(X,Y)}{{\sigma }_X{\sigma }_Y}$$
  为 $X$ 与 $Y$ 的 (线性) 相关系数，记为${\rho }_{xy}或Corr(X,Y)$。

相关系数${\rho }_{x,y}$的性质：

• 1. $-1⩽\mathrm{Corr}(X,Y)⩽1$， 或 $|\mathrm{Corr}(X,Y)|⩽1$。
• 2. $\mathrm{Corr}(X,Y)=\pm 1$ 的充要条件是 $X$ 与 $Y$ 间几乎处处有线性关系, 即存 在 $a(\ne 0)$ 与 $b$， 使得$$P(Y=aX+b)=1$$
• 3. 相关系数 $\mathrm{Corr}(X,Y)$ 刻画了 $X$ 与 $Y$ 之间的线性关系强弱， 因此也常称其为 “线性相关系数”。
• 4. 若 $\mathrm{Corr}(X,Y)=0$， 则称 $X$ 与 $Y$ 不相关。不相关是指 $X$ 与 $Y$ 之间没有线性关系， 但 $X$ 与 $Y$ 之间可能有其他的函数关系， 譬如平方关系、对数关系等。
• 5. 若 $\mathrm{Corr}(X,Y)=1$， 则称 $X$ 与 $Y$ 完全正相关； 若 $\mathrm{Corr}(X,Y)=-1$， 则称 $X$ 与 $Y$ 完全负相关。
• 6. 若 $0<|\mathrm{Corr}(X,Y)|<1$， 则称 $X$ 与 $Y$ 有 “一定程度” 的线性关系。 $|\mathrm{Corr}(X,Y)|$ 越接近于 1， 则线性相关程度越高； $|\mathrm{Corr}(X,Y)|$ 越接近于 0 ， 则线性相关程度越低。 而协方差看不出这一点， 若协方差很小， 而其两个标准差 ${\sigma }_X$ 和 ${\sigma }_Y$ 也很小， 则其比值就不一定很小。

  相关系数矩阵：类似于协方差矩阵，相关系数矩阵就是把协方差矩阵中每个元素替换成相关系数，具体来说就是：
$$\begin{array}{rl}& \mathrm{Corr}(X,Y)=\frac{\mathrm{Cov}(X,Y)}{\sqrt{\mathrm{Var}(X)}\sqrt{\mathrm{Var}(Y)}}=\frac{\mathrm{Cov}(X,Y)}{{\sigma }_X{\sigma }_Y}\\ =& \left(\begin{array}{cccc}1& \mathrm{Corr}\left(X_1,X_2\right)& \cdots & \mathrm{Corr}\left(X_1,X_n\right)\\ \mathrm{Corr}\left(X_2,X_1\right)& 1& \cdots & \mathrm{Corr}\left(X_2,X_n\right)\\ ⋮& ⋮& & ⋮\\ \mathrm{Corr}\left(X_n,X_1\right)& \mathrm{Corr}\left(X_n,X_2\right)& \cdots & 1\end{array}\right)\end{array}$$

python代码（求上一例题的相关系数矩阵）

from sympy import *
from sympy.abc import x, y 
f = (1 / 8) * (x + y)
fxy = Piecewise((f, And(x >= 0, x <= 2, y >= 0, y <= 2)), (0, True))
Ex = integrate(x * fxy, (x, -oo, oo), (y, -oo, oo))
Ey = integrate(y * fxy, (x, -oo, oo), (y, -oo, oo))
Exy = integrate(x * y * fxy, (x, -oo, oo), (y, -oo, oo))

cov_xy = Exy - Ex * Ey 
var_x = Ex_2 - (Ex ** 2)
var_y = Ey_2 - (Ey ** 2)
Matrix([[1, cov_xy/(var_x * var_y)], [cov_xy/(var_x * var_y), 1]])

$\left[\begin{array}{cc}1& -0.297520661157021\\ -0.297520661157021& 1\end{array}\right]$

1.7 随机变量序列的收敛状态：依概率收敛、依分布收敛

• 依概率收敛：设 { X n } \left\{X_{n}\right\} {Xn​} 为一随机变量序列， $X$ 为一随机变量， 如果对任意的 $\epsilon >0$， 有
  $$P\left(\left|X_n-X\right|⩾\epsilon \right)\to 0(n\to \infty )$$
  则称序列 { X n } \left\{X_{n}\right\} {Xn​} 依概率收敛于 $X$， 记作 $X_n\stackrel{P}{⟶}X$。

依概率收玫的含义是： $X_n$ 对 $X$ 的绝对偏差不小于任一给定量的可能性将随着 $n$增大而愈来愈小。或者说， 绝对偏差 $\left|X_n-X\right|$ 小于任一给定量的可能性将随着 $n$ 增大而愈来愈接近于 1 , 即$P\left(\left|X_n-X\right|⩾\epsilon \right)\to 0(n\to \infty )$等价于
$$P\left(\left|X_n-X\right|<\epsilon \right)\to 1(n\to \infty )$$
特别当 $X$ 为退化分布时， 即 $P(X=c)=1$（像概率p就是一个案例，频率不断趋近于一个常数p，这个p就是概率）， 则称序列 { X n } \left\{X_{n}\right\} {Xn​} 依概率收敛于 $c$， 即 $X_n\stackrel{P}{⟶}c$。

依概率收敛性质：设$X_n\stackrel{P}{⟶}a,Y_n\stackrel{P}{⟶}b$，又设函数$g(x,y)$在点$(a,b)$连续，则$$g(X_n,Y_n)\stackrel{P}{⟶}g(a,b)$$

• 依分布收敛：设随机变量 $X,X_1,X_2,\cdots$ 的分布函数分别为 $F(x),F_1(x),F_2(x),\cdots$。 若对 $F(x)$ 的任一连续点 $x$， 都有
  $$\underset{n\to \infty }{\lim }{F}_n(x)=F(x)$$
  则称 { F n ( x ) } \left\{F_{n}(x)\right\} {Fn​(x)} 弱收敛于 $F(x)$， 记作
  $$F_n(x)\stackrel{W}{⟶}F(x)$$
  也称相应的随机变量序列 { X n } \left\{X_{n}\right\} {Xn​} 按分布收敛于 $X$， 记作
  $$X_n\stackrel{L}{⟶}X$$

注：依概率收敛可推出依分布收敛，反之不成立。

1.8 大数定律

大数定律是叙述随机变量序列的前一些项的算数平均值在某种条件下收敛到这些项的均值的算数平均值。

  （弱大数）辛钦大数定理：设$X_1,X_2,\cdots$是相互独立，服从同一分布的随机变量序列，且$X_i,i=1,2,\cdots$的数学期望存在，作前$n$个变量的算数平均$\overline{X}=\frac{1}{n}{\sum }_{k=1}^n{X}_k$, 则对于任意的$\epsilon >0$，有
lim ⁡ n → ∞ P { ∣ 1 n ∑ k = 1 n X k − 1 n ∑ k = 1 n E ( X k ) ∣ < ε } = lim ⁡ n → ∞ P { ∣ X ‾ − E ( X ) ‾ ∣ < ε } = 1 \lim_{n \rightarrow \infty} P\{\left| \frac{1}{n} \sum_{k =1}^{n} X_{k} - \frac{1}{n} \sum_{k =1}^{n} E(X_{k}) \right| < \varepsilon\} = \lim_{n \rightarrow \infty} P\{\left| \overline{X} - \overline{E(X)}\right| < \varepsilon\}=1 n→∞lim​P{∣∣∣∣∣​n1​k=1∑n​Xk​−n1​k=1∑n​E(Xk​)∣∣∣∣∣​<ε}=n→∞lim​P{∣∣∣​X−E(X)​∣∣∣​<ε}=1
即 $\overline{X}\stackrel{P}{⟶}\overline{E(X)}$，若随机变量序列具有数学期望$E(X_k)=\mu (k=1,2,\cdots )$，则上式变为
lim ⁡ n → ∞ P { ∣ 1 n ∑ k = 1 n X k − μ ∣ < ε } = lim ⁡ n → ∞ P { ∣ X ‾ − μ ∣ < ε } = 1 \lim_{n \rightarrow \infty} P\{\left| \frac{1}{n} \sum_{k =1}^{n} X_{k} - \mu \right| < \varepsilon\} = \lim_{n \rightarrow \infty} P\{\left| \overline{X} - \mu \right| < \varepsilon\} = 1 n→∞lim​P{∣∣∣∣∣​n1​k=1∑n​Xk​−μ∣∣∣∣∣​<ε}=n→∞lim​P{∣∣​X−μ∣∣​<ε}=1
即 $\overline{X}\stackrel{P}{⟶}\mu$。

  伯努利大数定理：设$f_A$是$n$次独立重复实验中事件$A$发生的次数，$p$是事件$A$在每次试验中发生的概率，则对于任意正数$\epsilon >0$，有
lim ⁡ n → ∞ P { ∣ f A n − p ∣ < ε } = 1 \lim_{n \rightarrow \infty} P\{\left| \frac{f_{A}}{n} - p \right | < \varepsilon \} = 1 n→∞lim​P{∣∣∣∣​nfA​​−p∣∣∣∣​<ε}=1
或
lim ⁡ n → ∞ P { ∣ f A n − p ∣ ≥ ε } = 0 \lim_{n \rightarrow \infty} P\{\left| \frac{f_{A}}{n} - p \right | \ge \varepsilon \} = 0 n→∞lim​P{∣∣∣∣​nfA​​−p∣∣∣∣​≥ε}=0
即 $\frac{f_A}{n}\stackrel{P}{⟶}p$。

💡大数定理的条件：独立重复事件、重复次数足够多。

🔥例：使用蒙特卡洛模拟法求定积分

设 $0⩽f(x)⩽1$, 求 $f(x)$ 在 区间 $[0,1]$ 上的积分值
$$J={\int }_0^{1}f(x)\mathrm{d}x$$
🦊解:在正方形 { 0 ⩽ x ⩽ 1 , 0 ⩽ y ⩽ 1 } \{0 \leqslant x \leqslant 1,0 \leqslant y \leqslant 1\} {0⩽x⩽1,0⩽y⩽1}内均匀地投点$(x_i,y_i)$，投n个点，点越多越好。如果某个点$y_i\le f(x_i)$,则认为事件发生，我们计算满足$y_i\le f(x_i)$点的个数$S_n$，使用大数定律：频率稳定于概率，即：$\frac{S_n}{n}$就是积分值。

python代码（求解例题）

from scipy import integrate
from scipy.stats import uniform
import numpy as np
import matplotlib.pyplot as plt 
from random import choice

# 蒙特卡洛原理：将待积分函数看作是某个概率密度分布函数，根据分布函数的定义，那么对函数的积分就转化为求概率的过程。
x_arr = np.linspace(0, 1, 1000)
y_arr = x_arr ** 3

x_n = uniform.rvs(size=100)   
y_n = uniform.rvs(size=100)
plt.figure(figsize=(12, 6))
plt.stackplot(x_arr, y_arr, alpha=0.5, color='skyblue')
plt.scatter(x_n, y_n)
plt.text(1.0, 1.0, r'$y = x^3$')
plt.show()

class Montecarlo():
    def __init__(self, n):
        self.n = n    # 模拟的次数
        self.fn = 0 
    def solve(self, x, y, f):
        for i in range(self.n):
            a = choice(x)
            b = choice(y)
            if f(a) >= b:
                self.fn += 1
            else:
                continue
        res = self.fn / self.n
        print("蒙特卡洛方法(模拟{}次)计算的结果为：{}".format(self.n, res))

# 求解y = x^3 在区间[0, 1]上的积分值
x_1 = np.linspace(0, 1, 1000)
y_1 = np.linspace(0, 1, 1000)
def f(x):
    return x ** 3
int_f = integrate.quad(f, 0, 1)
print("函数y = x^3 的积分值为：{}".format(int_f[0]))
monte_10 = Montecarlo(10)
monte_10.solve(x_1, y_1, f)
monte_100 = Montecarlo(100)
monte_100.solve(x_1, y_1, f)
monte_1000 = Montecarlo(1000)
monte_1000.solve(x_1, y_1, f)
monte_10000 = Montecarlo(10000)
monte_10000.solve(x_1, y_1, f)

函数y = x^3 的积分值为：0.25
蒙特卡洛方法(模拟10次)计算的结果为：0.2
蒙特卡洛方法(模拟100次)计算的结果为：0.25
蒙特卡洛方法(模拟1000次)计算的结果为：0.247
蒙特卡洛方法(模拟10000次)计算的结果为：0.2521

1.9 中心极限定理

中心极限定理确定在什么条件下，大量随机变量之和的分布逼近于正态分布。

  独立同分布的中心极限定理：设随机变量$X_1,X_2,\cdots ,X_n,\cdots$相互独立，服从同一分布，且具有数学期望和方差：$E(X_k)=\mu ,D(X_k)={\sigma }^2(k=1,2,\cdots )$，则随机变量之和${\sum }_{k=1}^n{X}_k$的标准化变量
$$Y_n=\frac{{\sum }_{k=1}^n{X}_k-E({\sum }_{k=1}^n{X}_k)}{\sqrt{D({\sum }_{k=1}^n{X}_k)}}=\frac{{\sum }_{k=1}^n{X}_k-n\mu }{\sqrt{n}\sigma }$$
的分布函数$F_n(x)$对于任意$x$满足
$$\begin{array}{rl}\underset{n\to \infty }{\lim }{F}_n(X)& =\underset{n\to \infty }{\lim }P\left\{\frac{{\sum }_{k=1}^n{X}_k-n\mu }{\sqrt{n}\sigma }\le x\right\}\\ & ={\int }_{-\infty }^x\frac{1}{\sqrt{2\pi }}{e}^{-t^2/2}dt\\ & =\Phi (x)\end{array}$$
即 $Y_n\sim N(n\mu ,{\sigma }^2)$或$\overline{X}\sim N(\mu ,{\sigma }^2/n)$或$\frac{\overline{X}-\mu }{\sigma /\sqrt{n}}\sim N(0,1)$。

python代码（验证独立同分布的中心极限定理）

from sympy import * 
from scipy.stats import norm       # 正态分布
from scipy.stats import expon      # 指数分布
from scipy.stats import uniform     # 均匀分布
from scipy.stats import poisson     # 泊松分布
import scipy.stats as ss           # 0-1 分布(前文有重名的变量名，不便直接引入)  from scipy.stats import bernoulli
import matplotlib.pyplot as plt
import numpy as np
import random

# 1. n个服从正态分布随机变量的和 的分布

def norm_sum(mu, sigma, n, x_size):
    X = 0
    for i in range(n):
        x_i = norm(loc=mu[i], scale=sigma[i]).rvs(size=x_size)
        X += x_i
    plt.figure(figsize=(10, 8))    
    plt.hist(X, density=True, histtype='stepfilled', alpha=0.2)
    plt.title("Random variable number n={}".format(str(n)))
    plt.xlabel("x")
    plt.ylabel("p(x)")
    plt.show()

x_size = 10000
mu = np.zeros(x_size)
sigma = np.random.randint(0, 10, x_size)
norm_sum(mu, sigma, 10, x_size)
norm_sum(mu, sigma, 100, x_size)

# 2. n个服从均匀分布随机变量的和 的分布

def variable_sum(n, x_size, distribution):
    X = 0
    for i in range(n):
        x_i = distribution.rvs(x_size)
        X += x_i
    plt.figure(figsize=(10, 8))    
    plt.hist(X, density=True, histtype='stepfilled', alpha=0.2)
    plt.title("Random variable number n={}".format(str(n)))
    plt.xlabel("x")
    plt.ylabel("p(x)")
    plt.show()

x_size = 10000
loc = np.zeros(x_size)
scale = 5
distribution = uniform(loc=loc, scale=scale)
variable_sum(n=10, x_size=x_size, distribution=distribution)
variable_sum(n=100, x_size=x_size, distribution=distribution)

# 3. n个服从指数分布随机变量的和 的分布

x_size = 10000
lam = 5
distribution = expon(scale=1/lam)
variable_sum(n=10, x_size=x_size, distribution=distribution)
variable_sum(n=100, x_size=x_size, distribution=distribution)

# 4. n个服从泊松分布随机变量的和 的分布

x_size = 10000
mu = 5
distribution = poisson(mu=mu)
variable_sum(n=10, x_size=x_size, distribution=distribution)
variable_sum(n=100, x_size=x_size, distribution=distribution)

# 5. n个服从（0- 1）分布随机变量的和 的分布

def bernoulli_sum(n, x_size, p):
    X = 0
    for i in range(n):
        x_i = ss.bernoulli.rvs(p, size=x_size)
        X += x_i
    plt.figure(figsize=(10, 8))    
    plt.hist(X, density=True, histtype='stepfilled', alpha=0.2)
    plt.title("Random variable number n={}".format(str(n)))
    plt.xlabel("x")
    plt.ylabel("p(x)")
    plt.show()

x_size = 10000
p = 0.4
bernoulli_sum(n=10, x_size=x_size, p=p)
bernoulli_sum(n=100, x_size=x_size, p=p)

1.10 数学建模案例分析：投资组合分析

• 📕任务：GitModel公司是一家专业的投资银行，志在帮助客户更好地管理资产。客户手头上有一笔100万的资金，希望将这笔钱投入股票市场进行投资理财，投资人看中了两个股票$A$、$B$，股票分析师通过对股票$A$、$B$的历史数据分析发现：股票$A$的平均收益率近似服从$N(0.1,0.01)$，股票B的平均收益率近似服从$N(0.3,0.04)$。现在客户希望通过分析得出投资股票$A$、$B$的最佳组合（在预期收益率确定情况下最小风险时，需要投资$A$、$B$的份额）。

• 🦊解：由题意知：
  股票$A,B$的收益$X,Y$服从如下分布
  $$X\sim N(0.1,0.01),Y\sim N(0.3,0.04)$$
  设投资股票$A$的份额为$p$, 则投资股票$B$的份额为$(1-p)$，其中$(0\le p\le 1)$，总收益为$Z$，则
  $$Z=pX+(1-p)Y$$
  均值体现股票的平均收益率高低，方差体现股票的风险大小。则平均收益率为
  $$E(Z)=E(pX+(1-p)Y)=pE(X)+(1-p)E(Y)$$
  风险大小(目标函数)为
  $$\begin{array}{rl}D(Z)& =D(pX+(1-p)Y)\\ & =D(pX)+D((1-p)Y)+2Cov(pX,(1-p)Y)\\ & =p^{2}D(X)+(1-p{)}^{2}D(Y)+2p(1-p)Cov(X,Y)\end{array}$$
  由题意知，需要找到一个合适的$p\ast$，使得风险$D(Z)$最小
  $$venture=\underset{p^{\ast }}{\mathrm{argmin}}\left[D(Z)\right]=\underset{p^{\ast }}{\mathrm{argmin}}\left[p^{2}D(X)+(1-p{)}^{2}D(Y)+2p(1-p)Cov(X,Y)\right]$$

python代码（求解上题）

from sympy import *
from scipy.stats import norm 
from sympy.abc import p

X = norm(loc=0.1, scale=0.1)
Y = norm(loc=0.3, scale=0.2)
Dx, Dy = X.var(), Y.var()

def compute_risk(Dx, Dy, rho):
    # print("当相关系数rho={:.1f}时：".format(rho))
    Cov_xy = rho * sqrt(Dx) * sqrt(Dy)
    # 目标函数
    venture = p ** 2 * Dx + (1-p) ** 2 * Dy + 2 * p * (1-p) * Cov_xy
    # 求驻点
    venture_dp = diff(venture, p)
    p0 = solve(venture_dp)
    # print("\t驻点p0={:.2f}".format(p0[0]))
    # 判断在该点是否取得极小值
    venture_dp_2 = diff(venture_dp, p)
    temp = venture_dp_2.evalf(subs={p:p0[0]})
    if temp < 0:
        print("\t该点不是函数的极小值点。")      
    venture_min = venture.evalf(subs={p:p0[0]})
    # print("\t风险系数波动（方差）为：{:.4f}".format(venture_min))
    # print("\t此时投资股票A的份额应为：10w * {}, 投资股票B的份额应为：10w *{}".format(p0[0], 1-p0[0]))
    return p0[0], venture_min

rhos = [-1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
ps = []
ventures = []
for rho in rhos:
    temp1, temp2 = compute_risk(Dx, Dy, rho)
    ps.append(temp1)
    ventures.append(temp2)

for idx, p in enumerate(ps):
    r = (idx - 10) * 0.1
    if (p > 1) or (p < 0):
        print("相关系数rho={:.2f}时，投资方案无实际意义。".format(r))
    else:
        print("当相关系数rho={:.2f}时：".format(r))
        print("\t风险系数波动（方差）为：{:.4f}".format(ventures[idx]))
        print("\t此时投资股票A的份额应为：10w * {}, 投资股票B的份额应为：10w *{}".format(p, 1-p))

当相关系数rho=-1.00时：
	风险系数波动（方差）为：0.0000
	此时投资股票A的份额应为：10w * 0.666666666666667, 投资股票B的份额应为：10w *0.333333333333333
当相关系数rho=-0.90时：
	风险系数波动（方差）为：0.0009
	此时投资股票A的份额应为：10w * 0.674418604651163, 投资股票B的份额应为：10w *0.325581395348837
当相关系数rho=-0.80时：
	风险系数波动（方差）为：0.0018
	此时投资股票A的份额应为：10w * 0.682926829268293, 投资股票B的份额应为：10w *0.317073170731707
当相关系数rho=-0.70时：
	风险系数波动（方差）为：0.0026
	此时投资股票A的份额应为：10w * 0.692307692307692, 投资股票B的份额应为：10w *0.307692307692308
当相关系数rho=-0.60时：
	风险系数波动（方差）为：0.0035
	此时投资股票A的份额应为：10w * 0.702702702702703, 投资股票B的份额应为：10w *0.297297297297297
当相关系数rho=-0.50时：
	风险系数波动（方差）为：0.0043
	此时投资股票A的份额应为：10w * 0.714285714285714, 投资股票B的份额应为：10w *0.285714285714286
当相关系数rho=-0.40时：
	风险系数波动（方差）为：0.0051
	此时投资股票A的份额应为：10w * 0.727272727272727, 投资股票B的份额应为：10w *0.272727272727273
当相关系数rho=-0.30时：
	风险系数波动（方差）为：0.0059
	此时投资股票A的份额应为：10w * 0.741935483870968, 投资股票B的份额应为：10w *0.258064516129032
当相关系数rho=-0.20时：
	风险系数波动（方差）为：0.0066
	此时投资股票A的份额应为：10w * 0.758620689655172, 投资股票B的份额应为：10w *0.241379310344828
当相关系数rho=-0.10时：
	风险系数波动（方差）为：0.0073
	此时投资股票A的份额应为：10w * 0.777777777777778, 投资股票B的份额应为：10w *0.222222222222222
当相关系数rho=0.00时：
	风险系数波动（方差）为：0.0080
	此时投资股票A的份额应为：10w * 0.800000000000000, 投资股票B的份额应为：10w *0.200000000000000
当相关系数rho=0.10时：
	风险系数波动（方差）为：0.0086
	此时投资股票A的份额应为：10w * 0.826086956521739, 投资股票B的份额应为：10w *0.173913043478261
当相关系数rho=0.20时：
	风险系数波动（方差）为：0.0091
	此时投资股票A的份额应为：10w * 0.857142857142857, 投资股票B的份额应为：10w *0.142857142857143
当相关系数rho=0.30时：
	风险系数波动（方差）为：0.0096
	此时投资股票A的份额应为：10w * 0.894736842105263, 投资股票B的份额应为：10w *0.105263157894737
当相关系数rho=0.40时：
	风险系数波动（方差）为：0.0099
	此时投资股票A的份额应为：10w * 0.941176470588235, 投资股票B的份额应为：10w *0.0588235294117647
当相关系数rho=0.50时：
	风险系数波动（方差）为：0.0100
	此时投资股票A的份额应为：10w * 1.00000000000000, 投资股票B的份额应为：10w *0
相关系数rho=0.60时，投资方案无实际意义。
相关系数rho=0.70时，投资方案无实际意义。
相关系数rho=0.80时，投资方案无实际意义。
相关系数rho=0.90时，投资方案无实际意义。
相关系数rho=1.00时，投资方案无实际意义。