Box-muller Gaussian random Number generator & plot his

Question

我正在尝试使用Box-Muller方程用python编写代码，但我不知道如何开始！

这是我试图解决的示例：

• 在 900 keV 处观察到的峰值显示 FWHM 为 2 keV。使用下面列出的高斯采样方法，生成对应于 900 keV 峰值的 15,000 个计数，并节省采样能量。

• 创建并绘制箱宽度为 0.2 keV 的直方图，并与具有相同峰面积的高斯函数进行比较。

• 使用数据分析软件，尝试对蒙特卡罗数据进行高斯拟合，并查看结果是否足够接近峰值模型。

Box-Muller Gaussian抽样方法：[ 注意，下面的两个采样变量y1，y2是针对单位高斯分布(即mu=0，segma=1(。

y1 = (-2 ln r1)^1/2  * cos(2pi*r2)
y2 = (-2 ln r1)^1/2 * sin(2pi*r2)
(r1, r2: random numbers)}

有什么建议吗？

*更新*

我收到错误消息：

g1 = BoxMuller(v( 名称错误：未定义名称"v">

使用的代码是：

import random    
import matplotlib.pyplot as plt    
import numpy as np
def BoxMuller():    
r1 = np.random.randn(15000)*10    
r2 = np.random.randn(15000)    
a = 2.0 * np.pi * r2        
v = np.sqrt( -2.0*np.log(1.0 - r1)) * np.sin(a)        
u = np.sqrt( -2.0*np.log(1.0 - r1)) * np.cos(a)
g1 = BoxMuller(v)
g2 = BoxMuller(u)
q = 900.0 + g1*2.0
k = 900.0 + g2*2.0
plt.hist(q, k)
plt.show()

Answer 1

好吧，这是一个简单的实现来启动和修补

import math
import random
def BoxMuller():
r1 = random.random()
r2 = random.random()
a  = 2.0 * math.pi * r1
v  = math.sqrt( -2.0*math.log(1.0 - r2))
return (v * math.sin(a), v * math.cos(a))

g1, g2 = BoxMuller()
q = 900.0 + g1*2.0
...

更新

显然，给出了 FWHM，而不是 std.dev。要获得西格玛，必须将 FWHM 除以2*sqrt(2*log(2))~ 2.355。所以采样代码应该是

FWHM = 2.0
q = 900.0 + g1 * FWHM/2.355