Сгенерируйте случайные величины из распределения вероятностей

Я извлек некоторые переменные из своего набора данных 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

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


python numpy random math normal-distribution
person Pleastry    schedule 21.03.2020    source источник

Ответы (2)


arrow_upward
2
arrow_downward

Похоже, вы хотите сгенерировать данные на основе 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 PDF.

Вот пример кода. Я добавил здесь дополнительный уровень сложности, вычислив интеграл целевой функции от 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:]
person Mad Physicist    schedule 21.03.2020
comment
В моем случае я ничего не знаю о PDF во второй таблице, но это то, что я ожидаю на основе моих исходных данных. Поэтому я заменил данные в вашем ответе распределениями первой таблицы. x=[1, 2, 3, 4, 7, 13] -- p=[0.88, 0.07, 0.02, 0.01, 0.01, 0.01] Также была проблема с кодом, сопоставляющим скаляры с векторами, поэтому я изменил следующие строки: a = curve_fit(np.vectorize(model), x, p, 0.01) -- samples = np.ceil(expon.rvs(scale=1/a[0], size = 2000)).astype(int) В альтернативной модели нет необходимости в np.vectorize(). - person Pleastry; 23.03.2020
comment
@Черепаха. Звучит неплохо. Я бы сказал, что np.vectorize — ужасное решение для всего в этом контексте. Лучше на самом деле исправить причину ошибки. В чем ошибка? - person Mad Physicist; 23.03.2020
comment
@Черепаха. Кроме того, не забудьте в конечном итоге выбрать (и проголосовать) ответ. Он удалит вопрос из очереди без ответа, как только вы решите свою проблему. - person Mad Physicist; 23.03.2020
comment
Сгенерированный вывод — это то, что я ожидал бы, если бы использовал функцию экспоненциального затухания, но есть ли способ внести больше изменчивости в вывод? В настоящее время выходное распределение выглядит примерно так: values=[1, 2, 3, 4, 5] и counts=[1781, 185, 30, 3, 1], но я хотел бы, чтобы это было растянуто на более широкий диапазон значений, скажем, [1, 2, .., 15], если есть функция, которая могла бы это сделать, сохраняя при этом распределение, подобное тому, что первый стол. - person Pleastry; 23.03.2020
comment
исходная ошибка была TypeError: only size-1 arrays can be converted to Python scalars - person Pleastry; 23.03.2020
comment
@Черепаха. Весь смысл PDF в том, что это вероятность случайного распределения. - person Mad Physicist; 23.03.2020
comment
Давайте продолжим обсуждение в чате. - person Pleastry; 23.03.2020

arrow_upward
2
arrow_downward

Если у вас есть экспоненциальное затухание, базовое дискретное распределение вероятностей представляет собой геометрическое распределение. (Это дискретный аналог непрерывного экспоненциального распределения.) Такое геометрическое распределение использует параметр 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]
person JohanC    schedule 21.03.2020