在 python 中实施特定分配

Question

我想 return 1<l<10 概率 1/(2^(l-1))

我应该怎么做而不是：

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

等等

谢谢

Answer 1

这会很有趣...我对这些东西有点生疏，所以一个好的数学家可以解决我的推理。

要从公式生成分布，您首先需要进行一些积分并计算指定区间的累积密度函数。特别是我们需要开始计算归一化常数。

此积分给出，"k"：

累积密度函数的"meaning"为"what's the probability to obtain a certain number that belong to the interval I need?"。这道题可以换个角度看："the probability to take a number that is below or equal to 10 must be 1"。这导致以下等式有助于找到参数 "C"。请注意，第一个 therm 是 k，第二个 therm 是 2^(1-x) 的一般积分，其中我将 x 替换为 10。

解决这个问题我们终于得到了 CDF（同样，找到它的方法可能更简单）：

此时我们需要反转X的CDF。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中已经有截断指数，看一下here。

只需设置适当的比例和位置，即可完成工作。代码

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)

在 python 中实施特定分配

Implementing specific distribution in python

python

random

probability