Генерация случайных величин из распределения вероятностей

Question

Я извлек некоторые переменные из моего набора данных python, и я хочу сгенерировать больший набор данных из имеющихся у меня дистрибутивов. Проблема в том, что я пытаюсь внести некоторую изменчивость в новый набор данных, сохраняя при этом похожее поведение. Это пример моих извлеченных данных, который состоит из 400 наблюдений:

Value    Observation Count     Ratio of Entries
1        352                    0.88
2        28                     0.07
3        8                      0.02
4        4                      0.01
7        4                      0.01
13       4                      0.01

Теперь я пытаюсь использовать эту информацию для создания аналогичного набора данных с 2000 наблюдениями. Мне известны функции numpy.random.choice и random.choice, но я не хочу использовать точно такие же дистрибутивы. Вместо этого я хотел бы генерировать случайные величины (столбец значений) на основе распределения, но с большей изменчивостью. Пример того, как я хочу, чтобы мой больший набор данных выглядел следующим образом:

Value         Observation Count        Ratio of Entries
1             1763                     0.8815
2             151                      0.0755
3             32                       0.0160
4             19                       0.0095
5             10                       0.0050
6             8                        0.0040
7             2                        0.0010
8             4                        0.0020
9             2                        0.0010
10            3                        0.0015
11            1                        0.0005
12            1                        0.0005
13            1                        0.0005
14            2                        0.0010
15            1                        0.0005

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

Answer 1

Если у вас экспоненциальный спад, базовое дискретное распределение вероятностей имеет вид Геометрия c Распределение . (Это дискретный аналог непрерывного экспоненциального распределения .) Такое распределение геометрии c использует параметр p с вероятностью успеха одного испытания (например, смещенной монетой). Распределение описывает количество испытаний, необходимых для достижения одного успеха.

Ожидаемое среднее значение распределения 1/p. Таким образом, мы можем вычислить среднее значение наблюдений, чтобы оценить p.

. Функция является частью scipy как scipy.stats.geom. Чтобы получить образец распределения, используйте geom.rvs(estimated_p, size=2000).

. Вот код, демонстрирующий подход:

from scipy.stats import geom
import matplotlib.pyplot as plt
import numpy as np
from collections import defaultdict

observation_index = [1, 2, 3, 4, 7, 13]
observation_count = [352, 28, 8, 4, 4, 4]

observed_mean = sum([i * c for i, c in zip(observation_index, observation_count)]) / sum(observation_count)

estimated_p = 1 / observed_mean
print('observed_mean:', observed_mean)
print('estimated p:', estimated_p)

generated_values = geom.rvs(estimated_p, size=2000)
generated_dict = defaultdict(int)
for v in generated_values:
    generated_dict[v] += 1
generated_index = sorted(list (generated_dict.keys()))
generated_count = [generated_dict [i] for i in  generated_index]
print(generated_index)
print(generated_count)

Вывод:

observed_mean: 1.32
estimated p: 0.7575757575757576
new random sample:
    [1, 2, 3, 4, 5, 7]
    [1516, 365, 86, 26, 6, 1]

Answer 2

Похоже, вы хотите сгенерировать данные на основе PDF, описанного во второй таблице. PDF - это что-то вроде

0 for x <= B
A*exp(-A*(x-B)) for x > B

A определяет ширину вашего распределения, которая всегда будет нормализована, чтобы иметь область 1. B - это горизонтальное смещение, которое в вашем случае равно нулю , Вы можете сделать это целочисленным распределением, связав с помощью ceil.

CDF нормализованной убывающей экспоненты 1 - exp(-A*(x-B)). Как правило, простой способ создать пользовательский дистрибутив - это сгенерировать одинаковые числа и отобразить их через CDF.

К счастью, вам не придется этого делать, поскольку scipy.stats.expon уже обеспечивает реализацию, которую вы ищете. Все, что вам нужно сделать, это вписаться в данные в вашем последнем столбце, чтобы получить A (B явно равен нулю). Вы можете легко сделать это с помощью curve_fit. Помните, что A отображается на 1.0/scale на языке scipy.

Вот пример кода. Я добавил дополнительный уровень сложности, рассчитав интеграл целевой функции от n-1 до n для целочисленных входов, принимая во внимание биннинг для вас при выполнении подбора.

import numpy as np
from scipy.optimize import curve_fit
from scipy.stats import expon

def model(x, a):
    return np.exp(-a * (x - 1)) - exp(-a * x)
    #Alternnative:
    # return -np.diff(np.exp(-a * np.concatenate(([x[0] - 1], x))))

x = np.arange(1, 16)
p = np.array([0.8815, 0.0755, ..., 0.0010, 0.0005])
a = curve_fit(model, x, p, 0.01)
samples = np.ceil(expon.rvs(scale=1/a, size=2000)).astype(int)
samples[samples == 0] = 1
data = np.bincount(samples)[1:]