Руководство с использованием простого кода Python.
В моей предыдущей истории (https://oscarnieves100.medium.com/simulating-correlated-random-variables-in-python-c3947f2dbb10) я обсуждал, как моделировать коррелированные случайные числа из нормального распределения вероятностей. Мы сосредоточились на нормальном распределении, потому что оно хорошо известно и существуют стандартные методы выборки.
Но что, если вы хотите создать собственное распределение вероятностей? Что, если бы я сказал вам, что вам не нужно ограничиваться каким-либо из известных распределений вероятностей? Там, где есть воля, есть и способ, и то, что я собираюсь показать вам, очень простое, но очень мощное.
Сначала установим некоторые правила построения распределения вероятностей p (x):
- p (x) ≥ 0, то есть: p (x) должно быть положительным или нулевым для всех значений x, где x - действительное число. p (x) никогда не должно быть отрицательным, независимо от обстоятельств (отрицательные вероятности не имеют смысла).
- Интеграл от p (x) по всем возможным x должен быть равен единице. Если p (x) дискретно, то по определению сумма всех p (x) для всех x должна быть 1. Мы определяем это математически как Σp (x) = 1. В общем, мы можем гарантировать, что распределение вероятностей всегда подчиняется это свойство следующим образом: сначала мы создаем «ненормализованное» распределение g (x) в соответствии с некоторыми правилами по нашему выбору (например, g (x) = x² или как угодно, если соблюдается Правило № 1). Затем мы нормализуем его, разделив на сумму всех его значений, а именно: p (x) = g (x) / Σg (x). По определению это гарантирует, что Σp (x) = Σg (x) / Σg (x) = 1 всегда.
Теперь, когда мы разобрались с этим, давайте поиграемся с простым распределением: g (x) = (x-x0) ², где x0 - некоторое значение, которое мы определяем, которое сдвигает параболу вправо. Будем предполагать, что g (x) существует на интервале [a, b] и равна нулю вне этих границ.
Общая идея такова: сначала мы вычисляем массив X, который простирается от a до b, в идеале с равномерным интервалом. Затем мы вычисляем p (X) = g (X) / Σg (X). Мы будем называть этот массив P. Затем мы вычисляем совокупную сумму P, а именно C. Это означает, что для каждого значения в массиве X, C становится больше, потому что мы добавляем предыдущее значение в совокупную сумму. По определению максимальное значение C = 1, потому что ΣP = 1, поэтому C увеличивается от 0 до 1 некоторым образом, продиктованным P. Во-вторых, мы моделируем число из равномерного распределения, а именно:
U ~ Uniform (0,1) на отрезке [0, 1].
Затем мы умножаем U на окончательное значение C, называем это max (C). В качестве альтернативы мы можем просто выбрать последнее значение в C как C [-1], поскольку оно будет совпадать с max (C). Идея состоит в том, что для каждого значения в C мы проверяем, соответствует ли термин
U * макс (C) ›C
Это должно вернуть логический массив, который говорит True (1) или False (0) в каждой позиции массива C. Например, предположим, что C имеет в нем 5 значений. Тогда возможный логический массив из вычисления «U * max (C)› C »мог бы быть:
[Верно, Верно, Ложно, Ложно, Верно] → [1, 1, 0, 0, 1]
Затем мы суммируем значения в этом логическом массиве, как
Σ (U * макс (C) ›C)
Например, Σ ([1, 1, 0, 0, 1]) = 3. Это будет номер индекса, который мы выбираем, из которого мы выбираем значение в X, по сути, X [3]. Так почему это работает? По сути, каждый раз, когда вы выбираете U, вы получаете другой порядковый номер из Σ (U * max (C) ›C), который указывает на другое число в X. Мы будем называть эти числа R, так что
R = X [Σ (U * max (C) ›C)]
Поскольку член Σ (U * max (C) ›C) продиктован формой P, случайные числа, которые мы отбираем, будут следовать распределению P.
Вот как мы бы сделали это в Python:
# Arbitrary probability distribution # # Generates random numbers according to an arbitrary probability distribution # p(x) defined by the user on a domain [a,b] # # Author: Oscar A. Nieves # Last updated: July 30, 2021 import matplotlib.pyplot as plt import numpy as np import math as mt plt.close('all') #np.random.seed(0) # Set seed # Inputs a = 10 b = 25 x0 = (a+b)/2 samples = 1000 X = np.linspace(a,b,1001) # Define p(x) def g(x_var): return (x_var - x0)**2 P0 = g(X) # Normalize p(x) Psum = sum(P0) P = P0/Psum # Cumulative sum p(x) C = np.cumsum(P) # Compute numbers R = np.linspace(a,b,samples) for n in range(samples): index0 = mt.ceil( sum(sum(C[-1]*np.random.rand(1,1) > C)) ) R[n] = X[index0] # Generate histogram bins = 50 plt.subplot(1,2,1) plt.plot(X,P,color='black') plt.xlim([a,b]) plt.xlabel('X', fontsize=16) plt.ylabel('p(x)', fontsize=16) plt.subplot(1,2,2) plt.hist(R,bins) plt.xlim([a,b]) plt.xlabel('R', fontsize=16) plt.ylabel('frequency', fontsize=16)
Запуск этого кода дает следующие графики:
Слева у нас есть нормализованное распределение вероятностей p (x) (или P), а справа - гистограмма со случайно выбранными числами R. Обратите внимание, что здесь я создаю 1000 независимых выборок для R внутри цикла for. Интересно, что мы можем выбрать любую функцию g (x) по нашему выбору.
Например, давайте вместо этого используем g (x) = (x-x0) ⁸, который должен выглядеть намного шире в центре:
Вот и все, теперь вы знаете, как создавать любые случайные числа, какие захотите. Повеселись!