从截断的 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屋！