Реализация конкретного дистрибутива в Python

Question

я хочу вернуть 1<l<10 с вероятностью 1/(2^(l-1))

как мне это сделать, а не:

    x = random()
    if x < 0.5:
       return 2

и так далее

спасибо

Answer 1

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

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

Этот интеграл дает для «k»:

«Значение»кумулятивная функция плотности: «какова вероятность получить определенное число, которое принадлежит нужному мне интервалу?». Этот вопрос можно рассматривать по-другому: «вероятность взять число, которое меньше или равно 10, должно быть 1». Это приводит к следующему уравнению, которое помогает найти параметр «С». Обратите внимание, что первый терм - это k, второй терм - это общий интеграл от 2 ^ (1-x), где я должен заменить x на 10.

Решая это, мынаконец, достигните CDF (опять же, возможно, что найти его проще):

На этом этапе нам нужно повернуть CDF вспятьдля X. X теперь наш генератор случайных чисел от 0 до 1. Формула:

В коде Python я попробовал следующее:

import numpy as np
import matplotlib.pyplot as plt

a=[ 1-   np.log2(1-(1-2**(-9))*np.random.rand()) for i in range(10000)]

plt.hist(a, normed=True)

Имеет ли это смысл?

Answer 2

Хотя ответ @Fabrizio, вероятно, верен, существует гораздо более простой способ выполнить работу - вы хотите, чтобы это было усеченным экспоненциально, потому что ваш PDF выглядит как

PDF (x) ~ 2 ^-x = e ^{-x log (2)} .

В SciPy уже есть усеченная экспонента, посмотрите здесь .

Просто установите правильный масштаб и местоположение, и работа сделана. Код

import numpy as np
from scipy.stats import truncexpon
import matplotlib.pyplot as plt

vmin = 1.0
vmax = 10.0
scale=1.0/np.log(2.0)

r = truncexpon.rvs(b=(vmax-vmin)/scale, loc=vmin, scale=scale, size=100000)

print(np.min(r))
print(np.max(r))

plt.hist(r, bins=[1,2,3,4,5,6,7,8,9,10], density=True)

Гистограмма

И если вам нужно выбрать только целочисленные значения, в Numpy есть хорошая вспомогательная функция:ну, код ниже, график очень похож

#%%
import numpy as np
import matplotlib.pyplot as plt

vmin = 1
vmax = 10

v = np.arange(vmin+1, vmax, dtype=np.int64)
p  = np.asarray([1.0/2**(l-1) for l in range(vmin+1, vmax)]) # probabilities
p /= np.sum(p) # normalization

r = np.random.choice(v, size=100000, replace=True, p=p)

print(np.min(r))
print(np.max(r))

plt.hist(r, bins=[1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5], density=True)

Answer 3

В одну сторону:

x = random()
l =  2
total += 1/2**(l-1)
while x < total: 
    l += 1
    total += 1/2**(l - 1)
return l

, другой

x = 1 - random()
raturn math.ceil(1 - math.log(x, 2))