【Python】Numpy中的Gumbel分布和Logistic分布
文章目录
- 极值分
- Gumbel
- Logistic分布
极值分
设 X 1 , X 2 … , X n X_1,X_2\dots,X_n X1,X2…,Xn为从总体 F F F中抽出的独立同分布样本，且
M = max ( X 1 , … , X n ) , m = min ( X 1 , … , X n ) M=\max(X_1,\dots,X_n), m=\min(X_1,\dots,X_n) M=max(X1,…,Xn),m=min(X1,…,Xn)
若存在 C n > 0 C_n>0 Cn>0和 D n D_n Dn，使得 C n M + D n C_nM+D_n CnM+Dn按分布收敛于 G ( x ) G(x) G(x)，则此 G ( x ) G(x) G(x)为极大值分布，同理可定义极小值分布。Fisher和Tippett证明了极值分布只有三种形式，分别是
| I型 | G 1 ( x ) = exp ( − e − x ) G_1(x)=\exp(-e^{-x}) G1(x)=exp(−e−x) | Gumbel分布 |
| II型 | G 2 ( x ) = exp ( − x − α ) , x > 0 , α > 0 G_2(x)=\exp(-x^{-\alpha}), x>0, \alpha>0 G2(x)=exp(−x−α),x>0,α>0 | Fréchet分布 |
| III型 | G 3 ( x ) = exp ( − ( − x ) α ) , x < 0 , α > 0 G_3(x)=\exp(-(-x)^\alpha), x<0, \alpha>0 G3(x)=exp(−(−x)α),x<0,α>0 | Weibull分布 |
上面的三个表达式均未考虑x的偏移，若是存在存在位置 μ \mu μ和尺度 λ \lambda λ的偏移，可将 x x x替换为 x − μ λ \frac{x-\mu}{\lambda} λx−μ
Gumbel
在Numpy中，Gumbel分布被实现为gumbel([loc, scale])，其中loc和scale分别对应
μ
,
λ
\mu, \lambda
μ,λ，其具体表达式为
p ( x ) = exp [ − z − e − z ] , z = x − μ λ p(x)=\exp[{-z-e^{-z}}], z=\frac{x-\mu}{\lambda} p(x)=exp[−z−e−z],z=λx−μ
当 μ = 0 , λ = 1 \mu=0, \lambda=1 μ=0,λ=1时，得到标准极小值分布，其期望值为
E ( z ) = ∫ z e z e − e z d z = ∫ z e − e z d e z = ∫ 0 ∞ ln z e − z d z = − γ E(z)=\int ze^ze^{-e^z}\text dz=\int ze^{-e^z}\text de^z=\int_0^\infty\ln z e^{-z}\text dz=-\gamma E(z)=∫zeze−ezdz=∫ze−ezdez=∫0∞lnze−zdz=−γ
其中 γ ≈ − 0.577 \gamma\approx-0.577 γ≈−0.577，为欧拉常数，接下来可以尝试计算一下
import numpy as np
from numpy.random import gumbel
xs = gumbel(loc=0, scale=1, size=20000)
print(np.mean(xs))
# 0.5836924011336291
通过调整\lambda可以调整Gumbel分布的形状
import matplotlib.pyplot as plt
for lam in [2,1,0.5]:
xs = gumbel(scale=lam, size=20000)
plt.hist(xs, bins=200, label=f"lambda={lam}", alpha=0.5)
plt.legend()
plt.show()
如图所示，
可见\lambda越大，则Gumbel分布越宽。
Logistic分布
如果两个随机变量 X 1 X_1 X1和 X 2 X_2 X2均服从Gumbel分布，那么 X 1 − X 2 X_1-X_2 X1−X2服从Logistic分布。
如果把 X 1 X_1 X1和 X 2 X_2 X2理解为某一随机样本的最大值分布和最小值分布，那么Logistic分布也可以解释为，从指数分布的总体中，抽取容量为n的随机样本，当n趋于无穷大时，样本极差所处的分布状态，其概率密度函数为
p ( x ) = ( x − μ ) / s s ( 1 + exp [ − ( x − μ ) / s ] ) 2 p(x)=\frac{(x-\mu)/s}{s(1+\exp[-(x-\mu)/s])^2} p(x)=s(1+exp[−(x−μ)/s])2(x−μ)/s
其分布图像为
loc, scale = 10, 1
s = np.random.logistic(loc, scale, 10000)
plt.hist(s, bins=50)