python - Box-muller 高斯随机数生成器并绘制他的图
问题描述
我正在尝试使用 Box-Muller 方程在 python 中编写代码,但我不知道如何开始!
这是我要解决的示例:
在 900 keV 处观察到的峰显示出 2 keV 的 FWHM。使用下面列出的高斯采样方法,生成对应于 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()
解决方案
好吧,这是开始和修改的简单实现
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
推荐阅读
- android - 来自 GetCurrentMillis 的 Android Studio 日期比较
- javascript - 无法在 HTML 中加载 XMLHttpRequest 视频
- java - Weblogic/Java 未在与 IIS 的相互 SSL 集成中发送客户端证书
- azure - 在 Azure 中发布的 webapi 属于 PaaS 和 IaaS?
- asp.net - 找到 2 个具有非唯一 ID #__VIEWSTATE 的元素
- apache-spark - Oozie:超时后终止工作
- sql - 选择多个子查询
- c# - WPF 数据绑定到来自后面代码的 TreeView
- android - 拆分 Android 项目?
- r - 从R中的名称中提取数量元素