数学与编程——统计与编程

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

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

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

定理 Box-Muller变换，如果随机变量U1,U2独立且U1,U2∼Uniform[0,1]，则

Z0=−2lnU1−−−−−−−√cos(2πU2)Z1=−2lnU2−−−−−−−√cos(2πU1)

则Z0,Z1独立且服从标准正太分布。

import numpy as npimport matplotlib.pyplot as pltU1 = 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()

使用蒙特卡洛的方法逼近π

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