Как генерировать случайные числа с предопределенным распределением вероятностей?

Question

Я хотел бы реализовать функцию в Python (используя numpy), которая принимает математическую функцию (например, p(x) = e^(-x), как показано ниже) в качестве входных данных и генерирует случайные числа, которые распределяются в соответствии с вероятностью этой математической функции распределение. И мне нужно построить их так, чтобы мы могли увидеть распределение.

Мне действительно нужна функция генератора случайных чисел для следующих двух математических функций в качестве входных данных, но если она может принимать другие функции, почему бы и нет:

1) p(x) = e^(-x)
2) g(x) = (1/sqrt(2*pi)) * e^(-(x^2)/2)

Кто-нибудь знает, как это можно сделать в python?

Answer 1

Для простых дистрибутивов, таких как те, которые вам нужны, или если у вас есть возможность легко инвертировать CDF в закрытой форме, вы можете найти множество сэмплеров в NumPy, как правильно указано в ответе Оливье.

Для произвольных распределений вы можете использовать методы выборки Маркова-Цепи Монте-Карло.

Самым простым и, возможно, более простым для понимания вариантом этих алгоритмов является Метрополис выборка.

Основная идея выглядит так:

• начать со случайной точки x и сделать случайный шаг xnew = x + delta
• оценить желаемое распределение вероятностей в начальной точке p(x) и в новой p(xnew)
• если новая точка более вероятна p(xnew)/p(x) >= 1 принять ход
• если новая точка менее вероятна, случайным образом решите, принимать или отклонять, в зависимости от вероятности ¹ новая точка равна
• новый шаг с этой точки и повторите цикл

Это может быть показано, см., Например, Сокал ² , точки, выбранные этим методом, соответствуют распределению вероятности принятия.

Обширную реализацию методов Монте-Карло в Python можно найти в пакете PyMC3.

Пример реализации

Вот игрушечный пример, чтобы показать вам основную идею, никоим образом не подразумеваемую в качестве эталонной реализации. Пожалуйста, обратитесь к зрелым пакетам для любой серьезной работы.

def uniform_proposal(x, delta=2.0):
    return np.random.uniform(x - delta, x + delta)

def metropolis_sampler(p, nsamples, proposal=uniform_proposal):
    x = 1 # start somewhere

    for i in range(nsamples):
        trial = proposal(x) # random neighbour from the proposal distribution
        acceptance = p(trial)/p(x)

        # accept the move conditionally
        if np.random.uniform() < acceptance:
            x = trial

        yield x

Давайте посмотрим, работает ли он с некоторыми простыми дистрибутивами

Гауссова смесь

def gaussian(x, mu, sigma):
    return 1./sigma/np.sqrt(2*np.pi)*np.exp(-((x-mu)**2)/2./sigma/sigma)

p = lambda x: gaussian(x, 1, 0.3) + gaussian(x, -1, 0.1) + gaussian(x, 3, 0.2)
samples = list(metropolis_sampler(p, 100000))

Коши

def cauchy(x, mu, gamma):
    return 1./(np.pi*gamma*(1.+((x-mu)/gamma)**2))

p = lambda x: cauchy(x, -2, 0.5)
samples = list(metropolis_sampler(p, 100000))

Произвольные функции

Вам на самом деле не нужно выбирать из правильных распределений вероятностей. Возможно, вам просто придется ввести ограниченный домен, в котором вы можете выбрать случайные шаги ³

p = lambda x: np.sqrt(x)
samples = list(metropolis_sampler(p, 100000, domain=(0, 10)))

p = lambda x: (np.sin(x)/x)**2
samples = list(metropolis_sampler(p, 100000, domain=(-4*np.pi, 4*np.pi)))

Выводы

Еще слишком много можно сказать о распределении предложений, конвергенции, корреляции, эффективности, приложениях, байесовском формализме, других образцах MCMC и т. Д. Я не думаю, что это правильное место, и есть много гораздо лучшего материала, чем то, что я мог бы написать здесь, доступный в Интернете.

---

1. Идея в том, чтобы поддержать разведку, где вероятность выше, но все же смотреть на регионы с низкой вероятностью, поскольку они могут привести к другим пикам. Основополагающим является выбор дистрибутива Предложение , то есть, как вы выбираете новые точки для изучения. Слишком маленькие шаги могут ограничить вас в ограниченной области вашего распространения, слишком большие могут привести к очень неэффективному исследованию.

2. Физика ориентирована. Байесовский формализм (Метрополис-Гастингс) предпочтителен в наши дни, но ИМХО его немного сложнее понять новичкам. В Интернете доступно множество учебных пособий, например, см. это из Университета Дьюка.

3. Реализация, показанная не для того, чтобы не добавлять слишком много путаницы, но просто, вам просто нужно обернуть пробные шаги на краях домена или заставить желаемую функцию перейти в ноль за пределами домена.

Answer 2

NumPy предлагает широкий диапазон распределений вероятностей .

Первая функция - экспоненциальное распределение с параметром 1.

np.random.exponential(1)

Вторым является нормальное распределение со средним 0 и дисперсией 1.

np.random.normal(0, 1)

Обратите внимание, что в обоих случаях аргументы являются необязательными, поскольку они являются значениями по умолчанию для этих распределений.

В качестве идентификатора вы также можете найти эти дистрибутивы в модуле random как random.expovariate и random.gauss соответственно.

Более общие дистрибутивы

Хотя NumPy, вероятно, удовлетворит все ваши потребности, помните, что вы всегда можете вычислить обратную совокупную функцию распределения вашего распределения и входные значения из равномерного распределения .

inverse_cdf(np.random.uniform())

Например, если NumPy не предоставил экспоненциального распределения , вы можете сделать это.

exponential = lambda: -np.log(-np.random.uniform())

Если вы сталкиваетесь с дистрибутивами, которые CDF нелегко вычислить, рассмотрите отличный ответ filippo .