从截断的MaxwellBoltzmann分布中获取随机数

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的逆函数没有解析形式，我需要对我画的每个数字进行二等分。如果有一种方便的方法可以让我从截断的麦克斯韦-玻尔兹曼分布中以相当快的速度生成随机数，那就太好了

Answer 1

事实证明，有一种方法可以利用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的分析形式进行比较，两者非常匹配

Answer 2

在这里你可以做几件事

• 对于固定参数v_esc和v_0，n_0是一个常数，因此不需要在pdf方法中计算它
• 如果仅为SciPyrv_continuous子类定义PDF，则该类的rvs、mean等将非常缓慢，可能是因为该方法每次需要生成随机变量或计算统计数据时都需要集成PDF。如果速度很高，那么您将需要向maxwell_boltzmann_pdf添加一个使用自己的采样器的_rvs方法。（另请参见this question。）一种可能的方法是拒绝采样方法：在框中画一个数字，直到该框位于PDF中。它适用于具有有限域的任何有界PDF，只要您知道域和界是什么（界是域中f的最大值）。参见this question获取Python代码示例
• 如果您知道发行版的CDF，那么还有一些额外的技巧。其中之一是对连续分布进行抽样的相对较新的k-vector sampling方法。有两个阶段：设置阶段和采样阶段。设置阶段涉及通过根查找近似CDF的逆，采样阶段使用此近似生成随机数，以非常快的方式遵循分布，而无需进一步评估CDF。对于像这样的固定分布，如果您向我展示CDF，我可以预先计算使用该数据对分布进行采样所需的必要数据和代码。本质上，k-vector采样唯一重要的部分是根查找步骤
• 有关从任意分布进行采样的更多信息，请参见mysampling methods page