从截断的 Maxwell-Boltzmann 分布中获取随机数 [英] Getting random numbers from a truncated Maxwell-Boltzmann distribution
问题描述
我想使用截断的 Maxwell-Boltzmann 分布生成随机数.我知道 scipy 有内置的 Maxwell 随机变量,但没有它的截断版本(我也知道截断的正态分布,这在这里无关紧要).我尝试使用 rvs_continuous 编写自己的随机变量:
I would like to generate random numbers using a truncated Maxwell-Boltzmann distribution. I know that scipy has built-in Maxwell random variables, but there is no truncated version of it (I am also aware of a truncated normal distribution, which is irrelevant here). I have tried to write my own random variables using rvs_continuous:
import scipy.stats as st
class maxwell_boltzmann_pdf(st.rv_continuous):
def _pdf(self,x):
n_0 = np.power(np.pi,3/2)*np.square(v_0)*(v_0*erf(v_esc/v_0)-(2/np.sqrt(np.pi))*v_esc*np.exp(-np.square(v_esc/v_0)))
return (1/n_0)*(4*np.pi*np.square(x))*np.exp(-np.square(x/v_0))*np.heaviside(v_esc-x,0)
maxwell_boltzmann_cv = maxwell_boltzmann_pdf(a=0, b=v_esc, name='maxwell_boltzmann_pdf')
这正是我想要的,但对于我的目的来说太慢了(我正在做蒙特卡罗模拟),即使我在所有循环之外绘制了所有随机速度.我也想过使用逆变换采样方法,但是CDF的逆没有解析形式,我需要对我绘制的每个数字进行二等分.如果有一种方便的方法可以以不错的速度从截断的 Maxwell-Boltzmann 分布中生成随机数,那就太好了.
This does exactly what I want, but it is way too slow for my purpose (I am doing Monte Carlo simulations), even if I draw all the random velocities outside of all the loops. I have also thought of using Inverse transform sampling method, but the inverse of the CDF does not have an analytic form and I will need to do a bisection for every number I draw. It would be great if there is a convenient way for me to generate random numbers from a truncated Maxwell-Boltzmann distribution with decent speed.
推荐答案
事实证明,有一种方法可以使用 scipy 的 ppf 特征,通过逆变换采样方法生成截断的 Maxwell-Boltzmann 分布.我在这里发布代码以供将来参考.
It turns out that there is a way to generate a truncated Maxwell-Boltzmann distribution with the inverse transform sampling method using the ppf feature of scipy. I am posting the code here for future reference.
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import erf
from scipy.stats import maxwell
# parameters
v_0 = 220 #km/s
v_esc = 550 #km/s
N = 10000
# CDF(v_esc)
cdf_v_esc = maxwell.cdf(v_esc,scale=v_0/np.sqrt(2))
# pdf for the distribution
def f_MB(v_mag):
n_0 = np.power(np.pi,3/2)*np.square(v_0)*(v_0*erf(v_esc/v_0)-(2/np.sqrt(np.pi))*v_esc*np.exp(-np.square(v_esc/v_0)))
return (1/n_0)*(4*np.pi*np.square(v_mag))*np.exp(-np.square(v_mag/v_0))*np.heaviside(v_esc-v_mag,0)
# plot the pdf
x = np.arange(600)
y = [f_MB(i) for i in x]
plt.plot(x,y,label='pdf')
# use inverse transform sampling to get the truncated Maxwell-Boltzmann distribution
sample = maxwell.ppf(np.random.rand(N)*cdf_v_esc,scale=v_0/np.sqrt(2))
# plot the histogram of the samples
plt.hist(sample,bins=100,histtype='step',density=True,label='histogram')
plt.xlabel('v_mag')
plt.legend()
plt.show()
此代码生成所需的随机变量,并将其直方图与 pdf 的解析形式进行比较,两者非常匹配.
This code generates the required random variables and compare its histogram with the analytic form of the pdf, which matches with each other pretty well.
这篇关于从截断的 Maxwell-Boltzmann 分布中获取随机数的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!