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

逆变换采样

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

F−1(u)=inf{x|F(u)≥x}F−1(u)=inf{x|F(u)≥x} P(F−1(u)≤x)=P(u≤F(x))=P(u≤y)∣∣y=F(x)=y|y=F(x)=F(x)P(F−1(u)≤x)=P(u≤F(x))=P(u≤y)|y=F(x)=y|y=F(x)=F(x)

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

G(u)≈F(u)G(u)≈F(u)

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

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

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

GELU(x)=xΦ(x)=x2(1+erf(x√2))≈x2(1+tanh[√2π(x+0.044715x3)])GELU⁡(x)=xΦ(x)=x2(1+erf⁡(x2))≈x2(1+tanh⁡[2π(x+0.044715x3)]) Φ(x)≈12(1+tanh[√2π(x+0.044715x3)])Φ(x)≈12(1+tanh⁡[2π(x+0.044715x3)])

考虑到tanh(x)tanh⁡(x)的逆，

artanhx=12ln(1+x1−x)artanh⁡x=12ln⁡(1+x1−x) 0.044715x3+x−√π8ln(Φ(x)1−Φ(x))≈00.044715x3+x−π8ln⁡(Φ(x)1−Φ(x))≈0

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

0.044715x3+x−√π8ln(u1−u)=00.044715x3+x−π8ln⁡(u1−u)=0

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

x3+x0.044715−√π80.044715ln(u1−u)=0x3+x0.044715−π80.044715ln⁡(u1−u)=0

容易验证上式满足，

4p3+27q2>04p3+27q2>0

因此有实根，

x=3√−q2+√q24+p327+3√−q2−√q24+p327x=−q2+q24+p3273+−q2−q24+p3273

其中pp​是固定值，

p=10.044715≈22.3639p=10.044715≈22.3639

qq的取值和u∼U[0,1]u∼U[0,1]​相关，

q=−√π80.044715ln(u1−u)q=−π80.044715ln⁡(u1−u)

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

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

Φ(x)≈σ(1.702x)Φ(x)≈σ(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/