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

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

2 голосов
/ 24 апреля 2020

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

Ответы [ 2 ]

1 голос
/ 25 апреля 2020

Здесь можно сделать несколько вещей.

  • Для фиксированных параметров v_esc и v_0, n_0 - это константа, поэтому ее не нужно вычислять в pdf метод
  • Вы можете предварительно рассчитать значения PDF в различных точках от 0 до 550 (что является поддержкой распределения). Затем вы можете взять эти точки и их значения в качестве входных данных для кусочно-линейного распределения (которое похоже на то, которое найдено в C ++, но я не думаю, что это имеет SciPy), из которого вы затем берете выборку.
  • Если вы знаете CDF дистрибутива, то есть несколько дополнительных хитростей. Одним из них является относительно новый выборка k-вектора метод отбора непрерывного распределения. Есть две фазы: фаза настройки и фаза выборки. Этап настройки включает в себя аппроксимацию обратного значения CDF с помощью нахождения root, а на этапе выборки это приближение использует для генерации случайных чисел, которые очень быстро следуют за распределением, без дополнительной оценки CDF. Для такого фиксированного дистрибутива, как этот, если вы покажете мне CDF, я могу предварительно рассчитать необходимые данные и код, необходимый для выборки дистрибутива с использованием этих данных. По сути, единственная нетривиальная часть выборки k-вектора - это шаг root -обнаружения.
  • Более подробную информацию о выборке из произвольного распределения можно найти в моей выборке страница методов .
0 голосов
/ 25 апреля 2020

Оказывается, что существует способ генерирования усеченного распределения Максвелла-Больцмана с помощью метода выборки с обратным преобразованием с использованием функции 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 analyti c, которая очень хорошо соответствует друг другу.

...