从截断的Maxwell-Boltzmann分布中得到随机数

Question

我想用截断的Maxwell-Boltzmann分布来生成随机数。我知道has内置了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

结果表明，本文提出了一种利用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的解析形式进行比较，后者彼此匹配得很好。

Answer 2

这里有几件事你可以做。

• 对于固定参数v_esc和v_0，n_0是一个常数，因此不需要在pdf方法中进行计算。
• 如果只为SciPy rv_continuous子类定义PDF，则该类的rvs、mean等将非常慢，可能是因为该方法每次需要生成随机变量或计算统计数据时都需要集成PDF。如果速度很高，那么您就需要在maxwell_boltzmann_pdf中添加一个使用自己的采样器的_rvs方法。(另见这个问题。)一种可能的方法是拒绝抽样方法:在一个盒子里画一个数字，直到盒子落入PDF的范围。它适用于具有有限域的任何有界PDF，只要您知道域和界是什么(界是域中f的最大值)。例如，请参阅这个问题代码。
• 如果您知道发行版的CDF，那么还有一些额外的技巧。其中之一是用于连续分布抽样的相对新的https://arxiv.org/abs/2004.02339方法。有两个阶段:设置阶段和取样阶段。建立阶段包括通过根的查找来近似CDF的逆，并且采样阶段使用这种近似来生成随机数，这些随机数以非常快的方式跟随分布，而不必进一步地评估CDF。对于这样一个固定的发行版，如果您向我展示了CDF，我可以预先计算必要的数据和代码，使用这些数据对分布进行采样。从本质上讲，k向量采样的唯一重要部分是寻根步骤.
• 关于来自任意分布的抽样的更多信息，请在我的抽样方法页上找到。