随机模拟中，生成随机数是必须的步骤，我们之前的学习中了解了如何生成等概率的整数、如何生成[0,1)区间上均匀分布的随机数。我们也知道正态分布是现实中最常用的概率分布，那么，如何生成来自正态分布的随机数呢？很多数学家给出了近似的解决方案，比如著名的Box-Muller算法，计算方法是
𝑤=𝑠𝑖𝑛(2𝜋𝑣)(−2ln𝑢)1/2
其中𝑣和𝑢都是浮点型数字，是来自[0,1)区间上均匀分布的随机数。请书写一个函数rv()用于生成一个标准正态分布的随机数，这需要导入math和random模块。完成函数后，请将该函数加入到以下gaussian.py文件中，以导入模块的方法再次应用rv()函数
以下是gaussian.py文件内容：
import math
def pdf(x, mu=0.0, sigma=1.0):
"""正态分布概率密度函数，返回浮点型数字。正态分布密度函数表达式为
\phi(z) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-z^2/2}，其中z = (x - \mu) / \sigma
Keyword arguments:
x -- 浮点型数字，随机变量x的取值，
mu -- 随机变量x的均值，默认为0.0
sigma -- 随机变量x的标准差，默认为1.0
Return values:
函数返回一个浮点型数字。
"""
z = (x - mu) / sigma
density = 1 / (math.sqrt(2 * math.pi) * sigma) * math.exp(- z ** 2 / 2)
return density
def cdf(x, mu=0.0, sigma=1.0):
"""正态分布累积分布函数，返回浮点型数字。我们将使用泰勒展式近似计算累积分布函数的取值
\Phi(z) = \frac{1}{2} + \phi(z) \left( z + \frac{z^3}{3} +
\frac{z^5}{3\cdot 5} + \frac{z^7}{3\cdot 5 \cdot 7} + \cdots \right)
Keyword arguments:
x -- 浮点型数字，随机变量x的取值，
mu -- 随机变量x的均值，默认为0.0
sigma -- 随机变量x的标准差，默认为1.0
Return values:
函数返回一个浮点型数字。
"""
z = (x - mu) / sigma
if z < -8.0: return 0
if z > 8.0: return 1
total = 0.0
term = z
i = 3
epsilon = 1.0e-15
while abs(term) > epsilon:
total += term
term *= z ** 2 / i
i += 2
return 0.5 + total * pdf(z)
def reverse_cdf(p0, mu=0.0, sigma=1.0):
"""根据二分搜索计算正态分布给定概率p0对应的分位数
Keyword arguments:
p0 -- 概率取值，浮点型数字
mu -- 正态分布的均值，缺省值为0.0
sigma --正态分布标准差，缺省值为1.0
Return values:
返回正态分布的分位数，浮点型数字
"""
lower = -8.0 # 初始化区间下限
upper = 8.0 # 初始化区间上限
guess = (lower + upper) / 2.0
epsilon = 1e-10
while True:
# 只要猜测数值对应的概率与p0的差异大于epsilon就执行循环
prob = cdf(guess, mu=mu, sigma=1.0)
if abs(prob - p0) < epsilon:
break
return guess
elif prob > p0:
upper = guess
else :
lower = guess
guess = (lower + upper) / 2.0
return guess
def rv():
"""生成一个来自标准正态分布的随机数，
v, u 是来自0-1区间上的随机数，浮点型
返回一个标准正态分布的随机数，浮点型
"""
u = random.random()
v = random.random()
w = math.sin(2 * math.pi * v) * math.sqrt(-2 * math.log(u))
return w
def main():
x = 1.65
mu = 0.0
sigma = 1.0
print(cdf(x, mu, sigma))
print(rv())
if __name__ == '__main__':
main()
导入math和random模块，完成函数后，请将该函数加入到以下gaussian.py文件中，以导入模块的方法再次应用rv()函数
以下是修改后的gaussian.py文件内容：
import math
import random
def pdf(x, mu=0.0, sigma=1.0):
z = (x - mu) / sigma
density = 1 / (math.sqrt(2 * math.pi) * sigma) * math.exp(- z ** 2 / 2)
return density
def cdf(x, mu=0.0, sigma=1.0):
z = (x - mu) / sigma
if z < -8.0: return 0
if z > 8.0: return 1
total = 0.0
term = z
i = 3
epsilon = 1.0e-15
while abs(term) > epsilon:
total += term
term *= z ** 2 / i
i += 2
return 0.5 + total * pdf(z)
def reverse_cdf(p0, mu=0.0, sigma=1.0):
lower = -8.0 # 初始化区间下限
upper = 8.0 # 初始化区间上限
guess = (lower + upper) / 2.0
epsilon = 1e-10
while True:
# 只要猜测数值对应的概率与p0的差异大于epsilon就执行循环
prob = cdf(guess, mu=mu, sigma=1.0)
if abs(prob - p0) < epsilon:
break
return guess
elif prob > p0:
upper = guess
else :
lower = guess
guess = (lower + upper) / 2.0
return guess
def rv():
u = random.random()
v = random.random()
w = math.sin(2 * math.pi * v) * math.sqrt(-2 * math.log(u))
return w
if __name__ == '__main__':
main()
现在可以通过导入gaussian模块并调用rv()函数来生成标准正态分布的随机数了。
不知道你这个问题是否已经解决, 如果还没有解决的话:首先，Box-Muller算法是用来从标准正态分布中生成随机变量的方法。Box-Muller算法是一个典型的逆变换法，它将一个随机变量转换成一个服从正态分布的随机变量。其基本思想是将均匀分布的随机数转化成标准正态分布的随机数。具体步骤如下：
1.生成两个[0,1)区间上均匀分布的随机数 U1 和 U2
2.构造两个随机变量Z1和Z2：Z1 = sqrt(-2 log U1) cos(2piU2)，Z2 = sqrt(-2 log U1) sin(2piU2)
3.将Z1和Z2作为输出值，它们就是独立、标准正态分布的随机变量
针对这个问题，我们可以编写一个名为rv()的函数，该函数能够使用Box-Muller算法生成标准正态分布的随机数。
代码实现如下：
import numpy as np
def rv():
U1 = np.random.rand() # 生成[0,1)上的均匀分布随机数
U2 = np.random.rand() # 生成[0,1)上的均匀分布随机数
Z1 = np.sqrt(-2*np.log(U1))*np.cos(2*np.pi*U2) # 计算Z1
Z2 = np.sqrt(-2*np.log(U1))*np.sin(2*np.pi*U2) # 计算Z2
return Z1 # 返回Z1作为随机数
以上代码可以实现生成标准正态分布的随机数的功能，但是需要使用Python自带的随机数生成器。如果需要使用Box-Muller算法自己编写随机数生成器，则需要使用以下代码：
import numpy as np
def rv():
while True:
U1 = np.random.rand() # 生成[0,1)上的均匀分布随机数
U2 = np.random.rand() # 生成[0,1)上的均匀分布随机数
R2 = U1**2 + U2**2 # 计算半径R
if R2 < 1: # 判断半径是否在[0,1)之间
break
Z1 = np.sqrt(-2*np.log(R2)/R2)*U1 # 计算Z1
Z2 = np.sqrt(-2*np.log(R2)/R2)*U2 # 计算Z2
return Z1 # 返回Z1作为随机数
上面的代码中，我们使用了一个while循环来非常方便地实现了随机数生成，其中半径R的计算在本质上是一个拒绝采样过程。当半径R不在[0,1)之间时，我们就不断地生成新的随机数，直到得到一个满足条件的随机数为止。
最后，我们把rv()函数加入到gaussian.py文件中，以便在其他地方方便地调用。该文件的代码示例如下：
import numpy as np
def rv():
while True:
U1 = np.random.rand() # 生成[0,1)上的均匀分布随机数
U2 = np.random.rand() # 生成[0,1)上的均匀分布随机数
R2 = U1**2 + U2**2 # 计算半径R
if R2 < 1: # 判断半径是否在[0,1)之间
break
Z1 = np.sqrt(-2*np.log(R2)/R2)*U1 # 计算Z1
Z2 = np.sqrt(-2*np.log(R2)/R2)*U2 # 计算Z2
return Z1 # 返回Z1作为随机数
以上就是该问题的解决方案，我们通过实现Box-Muller算法，编写了rv()函数来生成标准正态分布的随机数，并且把这个函数加入到gaussian.py文件中，方便其他地方进行调用。