Руководство с использованием простого кода Python.

В моей предыдущей истории (https://oscarnieves100.medium.com/simulating-correlated-random-variables-in-python-c3947f2dbb10) я обсуждал, как моделировать коррелированные случайные числа из нормального распределения вероятностей. Мы сосредоточились на нормальном распределении, потому что оно хорошо известно и существуют стандартные методы выборки.

Но что, если вы хотите создать собственное распределение вероятностей? Что, если бы я сказал вам, что вам не нужно ограничиваться каким-либо из известных распределений вероятностей? Там, где есть воля, есть и способ, и то, что я собираюсь показать вам, очень простое, но очень мощное.

Сначала установим некоторые правила построения распределения вероятностей p (x):

  1. p (x) ≥ 0, то есть: p (x) должно быть положительным или нулевым для всех значений x, где x - действительное число. p (x) никогда не должно быть отрицательным, независимо от обстоятельств (отрицательные вероятности не имеют смысла).
  2. Интеграл от 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) ⁸, который должен выглядеть намного шире в центре:

Вот и все, теперь вы знаете, как создавать любые случайные числа, какие захотите. Повеселись!