数学与编程——统计与编程（均匀分布仿真高斯分布）

统计模拟中有一个重要的问题就是，给定一个概率分布$p(x)$，我们如何使用计算机生成（generate）服从该分布的样本。一般而言均匀分布$\textrm{Uniform}(0, 1)$的样本相对容易生成。

1. 使用均匀分布仿真高斯分布

均匀分布比较优良的性质是，我们常见的各种概率分布，无论是连续的还是离散的分布，都可以基于$\textrm{Uniform}(0, 1)$的样本生成，例如正太分布可以通过著名的$\textrm{Box-Muller}$变换得到，

定理 Box-Muller变换，如果随机变量$U_1,\,U_2$独立且$U_1,\,U_2\sim\textrm{Uniform}[0, 1]$，则

$$Z_0=\sqrt{-2\ln U_1}\cos(2\pi U_2)\\Z_1=\sqrt{-2\ln U_2}\cos(2\pi U_1)$$
则 $Z_0,\,Z_1$独立且服从标准正太分布。

import numpy as np
import matplotlib.pyplot as plt

U1 = np.random.uniform(0, 1, 1000)
U2 = np.random.uniform(0, 1, 1000)

Z0 = np.sqrt(-2*np.log(U1))*np.cos(2*np.pi*U2)

cnts, bins, _ = plt.hist(Z0, 30, normed=True)
plt.plot(bins, 1./np.sqrt(2*np.pi)*np.exp(-bins**2/2), c='r', lw=2)
plt.show()

2. 使用蒙特卡洛的方法逼近$\pi$

import numpy as np
N = 10000
cnt = 0
for i in range(N):
    x, y = np.random.uniform(0, 1, 2)
    if (x-1/2)**2+(y-1/2)**2 < 1/4:
        cnt += 1
print('sample times: {}, pi approx: {}'.format(N, 4*cnt/N))