从截断的 Maxwell-Boltzmann 分布中获取随机数

Question

我想使用截断的 Maxwell-Boltzmann 分布生成随机数。我知道 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')

这正是我想要的，但对于我的目的来说太慢了（我正在做蒙特卡洛模拟），即使我在所有循环之外绘制了所有随机速度。我也想过使用逆变换采样方法，但是 CDF 的逆没有解析形式，我需要对我绘制的每个数字进行二等分。如果有一种方便的方法可以让我从截断的 Maxwell-Boltzmann 分布中以相当快的速度生成随机数，那就太好了。

Answer 1

您可以在这里做几件事。

• 对于固定参数v_esc和v_0，n_0是一个常数，所以不需要在pdf方法中计算。
• 如果你只为一个SciPy的rv_continuous子类定义一个PDF，那么这个类的rvs、mean等会很慢，大概是因为方法每次需要集成PDF生成随机变量或计算统计量。如果速度非常重要，那么您将需要向maxwell_boltzmann_pdf 添加一个使用自己的采样器的_rvs 方法。 （另请参阅this question。）一种可能的方法是拒绝采样方法：在方框中画一个数字，直到方框落在 PDF 中。它适用于任何具有有限域的有界 PDF，只要您知道域和边界是什么（边界是域中 f 的最大值）。有关 Python 代码的示例，请参见 this question。
• 如果您知道发行版的 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 的解析形式进行比较，两者非常匹配。