Получение случайных чисел из усеченного распределения Максвелла-Больцмана

Я хотел бы генерировать случайные числа, используя усеченное распределение Максвелла-Больцмана. Я знаю, что scipy имеет встроенные случайные величины Максвелла, но его усеченной версии нет (я также знаю об усеченном нормальном распределении, которое здесь не имеет значения). Я попытался написать свои собственные случайные переменные, используя 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 не имеет аналитической формы, и мне нужно будет делать деление пополам для каждого числа, которое я рисую. Было бы здорово, если бы был удобный для меня способ генерировать случайные числа из усеченного распределения Максвелла-Больцмана с приличной скоростью.


python scipy random montecarlo
person vlee    schedule 24.04.2020    source источник
comment
Каковы типичные диапазоны параметров?   -  person Peter O.    schedule 24.04.2020
comment
Я использую v_0 = 220 км/с и v_esc = 550 км/с.   -  person vlee    schedule 24.04.2020
comment
Это единственные значения, которые вы используете для v_0 и v_esc?   -  person Peter O.    schedule 24.04.2020
comment
Да (по крайней мере, для текущей симуляции, над которой я работаю)!   -  person vlee    schedule 24.04.2020

Ответы (2)


arrow_upward
1
arrow_downward

Есть несколько вещей, которые вы можете сделать здесь.

  • Для фиксированных параметров v_esc и v_0 n_0 является константой, поэтому ее не нужно вычислять в методе pdf.
  • Если вы определяете только PDF для подкласса SciPy rv_continuous, то классы rvs, mean и т. д. будут очень медленными, предположительно потому, что метод должен интегрировать PDF каждый раз, когда ему нужно создать случайную переменную или вычислить статистику. Если скорость имеет большое значение, вам необходимо добавить к maxwell_boltzmann_pdf метод _rvs, использующий собственный сэмплер. (См. также этот вопрос.) Одним из возможных методов является метод отклонение выборки: нарисуйте число в поле, пока поле не попадет в PDF. Он работает для любого ограниченного PDF-файла с конечным доменом, если вы знаете, что такое домен и граница (границей является максимальное значение f в домене). См. этот вопрос, например код Python.
  • Если вы знаете CDF дистрибутива, то есть некоторые дополнительные хитрости. Одним из них является относительно новый метод k-векторной выборки для выборки непрерывная раздача. Есть две фазы: фаза установки и фаза выборки. Фаза настройки включает в себя аппроксимацию обратной функции CDF посредством нахождения корня, а фаза выборки использует эту аппроксимацию для генерации случайных чисел, которые очень быстро следуют распределению без необходимости дальнейшей оценки функции CDF. Для такого фиксированного дистрибутива, как этот, если вы покажете мне CDF, я смогу предварительно рассчитать необходимые данные и код, необходимый для выборки дистрибутива с использованием этих данных. По сути, единственной нетривиальной частью выборки k-векторов является этап поиска корня.
  • Дополнительную информацию о выборке из произвольного дистрибутива можно найти на моей странице методов выборки.
person Peter O.    schedule 24.04.2020

arrow_upward
0
arrow_downward

Оказывается, есть способ создать усеченное распределение Максвелла-Больцмана с помощью метода выборки с обратным преобразованием, используя функцию ppf в scipy. Я отправляю код здесь для дальнейшего использования.

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, которые довольно хорошо совпадают друг с другом.

person vlee    schedule 25.04.2020