一种基于光滑逼近的正态分布采样法

一种基于光滑逼近正态分布的累积分布函数的采样法

逆变换采样

累积分布函数（CDF）是个单调函数，那么累积分布函数的反函数为，

$$F^{-1}(u) = \inf \left \{ x |F(u) \geq x \right \}$$

证明，

$$\begin{align} P(F^{-1}(u) \le x) &= P(u \le F(x)) \newline &= P(u \le y) \Big|_{y = F(x)} \newline &= y|_{y = F(x)} \newline &= F(x) \end{align}$$

这里有一个细节需要说明，由于$u$采样自均匀分布，因此容易知道$P(u \le y) = y$。这一点是以上推导的关键细节。对于很多分布来说，逆函数$F^{-1}(u)$并不容易计算，因此很多情况下无法直接使用逆变换采样。一个折中的思路是寻找$F(u)$的近似，

$$G(u) \approx F(u)$$

使得$G(u)$​容易求逆。于是$G^{-1}(u)$​​近似服从目标分布，$G(u)$逼近$F(u)$的误差越小，采样越接近目标分布。

利用这个思路，这里分享一种最近想到的正态分布采样方法。

原理

根据论文Gaussian Error Linear Units (GELUs)中的结论，

$$\begin{align} \operatorname{GELU}(x) &= x \Phi(x) \newline &= \frac{x}{2} \Big( 1 + \operatorname{erf}(\frac{x}{\sqrt{2}}) \Big) \newline &\approx \frac{x}{2} \Big( 1 + \tanh \left[ \sqrt{\frac{2}{\pi}} (x + 0.044715x^3)\right] \Big) \end{align}$$

于是有，

$$\Phi(x) \approx \frac{1}{2} \Big( 1 + \tanh \left[ \sqrt{\frac{2}{\pi}} (x + 0.044715x^3)\right] \Big)$$

考虑到$\tanh(x)$的逆，

$$\operatorname {artanh} x=\frac {1}{2}\ln \left({\frac {1+x}{1-x}}\right)$$

于是有，

$$0.044715x^3 + x - \sqrt{\frac{\pi}{8}}\ln \left({\frac {\Phi(x)}{1-\Phi(x) }}\right) \approx 0$$

根据逆变换采样原理，另$\Phi(x) = u \sim U[0, 1]$，求解一元三次方程，

$$0.044715x^3 + x - \sqrt{\frac{\pi}{8}}\ln \left({\frac {u}{1-u }}\right) = 0$$

获得的实根服从正态分布（近似）。接下来就是解一元三次方程的问题，化成$x^3 + px + q=0$​​的形式，于是有，

$$x^3 + \frac{x}{0.044715} - \frac{\sqrt{\frac{\pi}{8}}}{0.044715}\ln \left({\frac {u}{1-u }}\right) = 0$$

容易验证上式满足，

$$4p^{3}+27q^{2}>0$$

因此有实根，

$$x = \sqrt[{3}]{-\frac {q}{2}+\sqrt {\frac {q^{2}}{4}+\frac {p^{3}}{27}}}+\sqrt[{3}]{-\frac {q}{2}-\sqrt {\frac {q^{2}}{4}+\frac {p^{3}}{27}}}$$

其中$p$​是固定值，

$$p = \frac{1}{0.044715} \approx 22.3639$$

$q$的取值和$u \sim U[0, 1]$​相关，

$$q = - \frac{\sqrt{\frac{\pi}{8}}}{0.044715}\ln \left({\frac {u}{1-u }}\right)$$

求得的$x$​的取值​服从正态分布（近似）。

此外，论文Gaussian Error Linear Units (GELUs)中还有近似，

$$\Phi(x) \approx \sigma(1.702x)$$

不过该逼近的效果差于前者。

实现

Python实现如下（这里提供两种实现），

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

def normal_generator():
    # 实现一
    while True:
        u = np.random.uniform()
        x = np.sqrt(np.pi / 8) * np.log(u / (1 - u))
        r = np.roots([0.044715, 0, 1, -x])
        yield r[-1].real
        
def normal_generator2():
    # 实现二
    p = 1 / 0.044715
    t = -np.sqrt(np.pi / 8) / 0.044715
    while True:
        u = np.random.uniform()
        u = np.log(u / (1 - u))
        q = t * u
        d = np.sqrt(q ** 2 / 4 + p ** 3 / 27)
        r = np.cbrt(-q / 2 + d) + np.cbrt(-q / 2 - d)
        yield r

size = 10000
gen = normal_generator()
ns = np.array([next(gen) for _ in range(size)])
print(stats.normaltest(ns))
plt.hist(ns, bins=100, density=True)

mu = 0
sigma = 1
x = np.linspace(np.min(ns), np.max(ns), size)
y = stats.norm.pdf(x, mu, sigma)
plt.plot(x, y, 'r')
plt.show()

正态分布检验，可以看到pvalue还是比较显著的，正态效果不易于Numpy自带的正态分布采样方法。

NormaltestResult(statistic=0.9375257930807481, pvalue=0.6257759392078861)

注意，这里正太检验使用的方法是Pearson方法，具体就是偏度系数和峰度系数测试正态性。

可视化结果，

总结

本文分享一种基于逆变换采样与函数光滑近似的正态分布采样方法。

转载请包括本文地址：https://allenwind.github.io/blog/14990
更多文章请参考：https://allenwind.github.io/blog/archives/