Неправильное распределение из выборки обратного преобразования (CDF)

Question

Я пытаюсь смоделировать геометрическое распределение, используя метод Inverse CDF, однако я получаю немного неправильные результаты, и я не уверен, почему.

Точнее, геометрическое распределение с коэффициентом формыр = 0,8, должен иметь следующие характеристики:

mean: 1.25 
variance: 0.31

Однако, запустив приведенный ниже код, я получаю:

mean: 0.6224363901913519
var: 0.391813011265263
[Finished in 0.3s]

Как вы видите, я получаю дико другоесреднее значение по сравнению с ожидаемым.

np.log (iform [i]) / np.log (1-p) является результатом решения уравнения: F (X) = R дляX в терминах R, F (X) = CDF геометрического распределения = 1 - (1 - p) ^ k.

R - равномерное распределение по интервалу (0,1).

Таким образом, его решение приводит к следующему:

X = ln (1-R) ​​/ ln (1-p)

Однако, поскольку оба1-R и R равномерно распределены по (0,1), мы можем сделать следующее упрощение:

X = ln (R) / ln (1-p)

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

import numpy as np

n = 10000
p = 0.8
geo_dist = np.zeros(n,dtype = np.float64)
uniform = np.random.uniform(0, 1, n)
for i in range(n):
    geo_dist[i] = np.log(uniform[i])/np.log(1-p)
print("mean: " +str(geo_dist.mean()))
print("var: " +str(geo_dist.var())) 

Я попытался увеличить точность вычислений с помощью np.float64 в отчаянной попытке исправить то, что должно быть тривиальным скриптом, но безрезультатно.

Я также попытался сгенерировать равномерное распределение, используя scipyiform.rvs () вместо np.uniform, и проблема сохраняется.

Если p = 0,5:

expected mean: 2
expected variance : 2

код, который я написал, имеет следующий результат:

mean: 1.4440009653569306
var: 2.0421079966161093
[Finished in 0.3s]

У кого-нибудь есть идеи, почему это не работает?Спасибо.

Answer 1

Вы на самом деле производите непрерывную выборку экспоненциальное распределение с лямбдой, равной -1 / ln (1-p)

Хорошо, вот код с правильной выборкой, применяется потолокк экспоненциальному выводу

import numpy as np

N = 100000
p = 0.8

q = np.random.random(N)
g = np.ceil(np.log(1.0 - q)/np.log(1.0-p))

print(np.mean(g))
print(np.var(g))

, который печатает

1.25055
0.3146946975

Обратите внимание:

1. Вам лучше использовать возможности векторизации NumPy без явногопетли

2. Замена (1-R) -> R для R, отобранных из U (0,1), не работает для NumPy RNG - он возвращает значения в полузакрытом диапазоне [0 ...1), что означает, что вы можете время от времени получать исключения log (0) и FP.