从截断的麦克斯韦-玻尔兹曼分布中获取随机数

Question

我想使用截断的麦克斯韦-玻尔兹曼分布生成随机数。我知道 scipy 有内置的 Maxwell 随机变量，但没有它的截断版本（我也知道截断的正态分布，这里无关紧要）。我尝试使用 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')

这完全符合我的要求，但对于我的目的来说太慢了（我正在做 Monte Carlo 模拟），即使我在所有循环之外绘制了所有随机速度。我也想过使用Inverse transform sampling方法，但是CDF的逆没有解析形式，我需要对我绘制的每个数字进行二等分。如果有一种方便的方法可以让我以合适的速度从截断的麦克斯韦-玻尔兹曼分布中生成随机数，那就太好了。

Answer 1

您可以在这里做几件事。

对于固定参数v_esc和v_0，n_0是一个常数，所以不需要在pdf方法中计算。
如果您仅为 SciPy rv_continuous subclass 定义 PDF，则 class 的 rvs、mean ,等等会很慢，大概是因为该方法每次需要生成随机变量或计算统计量时都需要对PDF进行积分。如果速度非常重要，那么您需要向 maxwell_boltzmann_pdf 添加一个使用自己的采样器的 _rvs 方法。（另见 .) One possible method is the rejection sampling method: Draw a number in a box until the box falls within the PDF. It works for any bounded PDF with a finite domain, as long as you know what the domain and bound are (the bound is the maximum value of f in the domain). See this question 例如 Python 代码。
如果您知道发行版的 CDF，那么还有一些额外的技巧。其中之一是用于对连续分布进行采样的相对较新的 k-vector sampling 方法。有两个阶段：设置阶段和采样阶段。设置阶段涉及通过求根来近似 CDF 的逆，而采样阶段使用此近似值以非常快的方式生成服从分布的随机数，而无需进一步评估 CDF。对于像这样的固定分布，如果您向我展示 CDF，我可以预先计算必要的数据以及使用该数据对分布进行采样所需的代码。本质上，k 向量采样 唯一重要的部分是求根步骤。
有关从任意分布抽样的更多信息，请访问我的 sampling methods page。

Answer 2

原来有一种方法可以利用scipy的ppf特征，用逆变换采样的方法生成截断的Maxwell-Boltzmann分布。我在这里发布代码以供将来参考。

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 的解析形式进行比较，两者非常匹配。

从截断的麦克斯韦-玻尔兹曼分布中获取随机数

Getting random numbers from a truncated Maxwell-Boltzmann distribution

python

random

scipy

montecarlo