Box-muller 高斯随机数生成器并绘制他的图

Question

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

这是我要解决的示例：

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

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

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

Box-Muller 高斯采样方法： [注意，下面的两个采样变量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) NameError: name 'v' is not defined

使用的代码是：

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。要获得 sigma，必须将 FWHM 除以 2*sqrt(2*log(2)) ~ 2.355。所以采样代码应该是

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