概
感覺這個分布的含義很有用啊, 能預測‘最大', 也就是自然災害, 太牛了.
主要內容
定義
[Gumbel distribution-wiki](Gumbel distribution - Wikipedia)
其分布函數和概率密度函數分別為:
\[F(x; \mu, \beta) = e^{-e^{-(x-\mu)/\beta}}, \quad f(x;\mu,\beta) = \frac{1}{\beta} e^{-[e^{-(x-\mu)/\beta}+(x-\mu) / \beta]} \]
標准Gumbel分布(即\(\mu=0, \beta=1\)):
\[F(x) = e^{-e^{-x}}, \quad f(x) = e^{-(x+e^{-x})}. \]
從Gumbel分布中采樣, 只需:
\[x = F^{-1}(u) = \mu - \beta \ln (-\ln(u)), \quad u \sim \mathrm{Uniform}(0, 1). \]
proof:
\[P(F^{-1}(u) \le x) = P(u \le F(x)) = F(x), \]
故\(F^{-1}(u)\)的分布函數就是\(F(x)\).
\[\mathbb{E} [x] = \mu + \gamma \cdot \beta, \]
其中 \(\gamma\)是Euler-Mascherorni constant.
Gumbel-Max trick
假設我們有一個離散的分布\([\pi_1, \pi_2, \cdots, \pi_k]\)共\(k\)類, \(\pi_i\)表示為第\(i\)類的概率, 則從該分布中采樣\(z\)等價於
\[z = \arg \max_i [g_i + \log \pi_i], \quad g_i \sim \mathrm{Gumbel}(0, 1), \mathrm{i.i.d}. \]
proof:
\[P(z=i) = P(g_i + \log \pi_i \ge \max \{g_j + \log \pi_j\}_{j\not=i}) = \int_{-\infty}^{+\infty} p(x) P(x+\log \pi_i \ge \{g_j + \log \pi_j\}_{j\not=i}) \mathrm{d}x. \]
又
\[P(x+\log \pi_i \ge \{g_j + \log \pi_j\}_{j\not=i}) = \prod_{j\not=i} P(g_j \le x + \log\pi_i - \log \pi_j) = e^{-e^{-x} \cdot \frac{1 - \pi_i}{\pi_i}}, \]
帶入計算得:
\[\begin{array}{ll} P(z=i) & = \int_{-\infty}^{+\infty} e^{-(x+e^{-x} \cdot \frac{1}{\pi_i})} \mathrm{d}x \\ & = \int_{-\infty}^{+\infty} \pi_i \cdot e^{-[(x-\log\frac{1}{\pi_i})+e^{-(x - \log \frac{1}{\pi_i})}]} \mathrm{d}x \\ & = \pi_i. \end{array} \]
Gumbel trick 用於歸一化
我們時常會碰到這樣的問題:
\[p(x;\theta) = \frac{f(x;\theta)}{Z}, \]
其中\(Z=\sum_{i=1}^K f(x_i;\theta)\) 是歸一化常數, 那么怎么計算\(Z\)呢?
構建隨機變量\(T\):
\[T = \max_i [\ln f(x_i) + g_i], \quad g_i \sim \mathrm{Gumbel}(-c, 1), \mathrm{i.i.d.} \]
則
\[T \sim \mathrm{Gumbel}(-c + \ln Z) \]
proof:
\[P(T \le t) = P(\max_i [\ln f(x_i) + g_i] \le t) = \prod_{i} P(g_i \le t - \ln f(x_i)) = e^{-e^{-(t+c-\ln Z)}} = F(t;-c+\ln Z ,1). \]
因為
\[\mathbb{E}[T] = -c + \ln Z + \gamma, \]
故我們只需估計\(\mathbb{E}[T] \approx \sum_j T_j\) 即可估計\(Z\)
\[Z = \exp (\sum_{j}T_j + c - \gamma). \]
所以必須要求離散的\(x\)?
代碼
[scipy-gumbel](scipy.stats.gumbel_r — SciPy v1.6.3 Reference Guide)
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gumbel_r
fig, ax = plt.subplots(1, 1)
# mean, var, skew, kurt = gumbel_r.stats(moments='mvsk')
# print(mean, var, skew, kurt)
x = np.linspace(gumbel_r.ppf(0.01), gumbel_r.ppf(0.99), 100)
ax.plot(x, gumbel_r.pdf(x), 'r-', lw=5, alpha=0.6, label="gumbel_r pdf")
r = gumbel_r.rvs(size=1000, loc=0, scale=1)
ax.hist(r, density=True, histtype="stepfilled", alpha=0.2)
ax.legend(loc='best', frameon=False)
plt.show()
