在Python中如何对截断整数幂律进行抽样?

7

如果我想对截断整数幂律进行采样,我该使用Python中的哪个函数?

也就是说,给定两个参数am,在范围[1,m)内生成一个随机整数x,其遵循与1/x^a成比例的分布。

我已经在numpy.random中搜索过,但没有找到这个分布。


为什么不使用内置的幂律分布进行拒绝抽样呢? - Zach H
3个回答

5
据我所知,NumPy和Scipy都没有为您定义这个分布。但是,使用SciPy可以使用scipy.rv_discrete轻松地定义自己的离散分布函数:
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

def truncated_power_law(a, m):
    x = np.arange(1, m+1, dtype='float')
    pmf = 1/x**a
    pmf /= pmf.sum()
    return stats.rv_discrete(values=(range(1, m+1), pmf))

a, m = 2, 10
d = truncated_power_law(a=a, m=m)

N = 10**4
sample = d.rvs(size=N)

plt.hist(sample, bins=np.arange(m)+0.5)
plt.show()

enter image description here


看起来你正在将pmf集成,就好像它是连续的一样,并且在1和2之间取面积以得出p(1),在2和3之间取面积以得出p(2),等等,是这样吗?如果是这样,在你的例子中,我认为你需要模仿Spinal Tap并到11去获取p(10)。你的“const”将通过在分母中有“(m+1)**k”进行调整。或者我误解了? - pjs
@pjs:我将pdf视为连续函数1/x**a。因此,没有在区间[1,2]、[2,3]等上进行积分。但是,我手动积分找到了const_ppf的公式,即cdf的反函数。我认为我做对了,但我可能错了。(我尝试了您的建议,但它将域移位到[1,11],所以如果我理解正确,这不符合基本的合理性检查。)顺便问一下,Spinal Tap在这里指的是什么? - unutbu
我不是Pythonista,所以无法直接检查您的结果,但我已经针对a,m = 2,10进行了直接计算,p(1)应该为0.6452579827864142。这是您得到的结果吗? - pjs
@pjs:我手动计算得到了相同的结果(0.555...),通过对函数1/x**2x=1x=2进行积分,并进行归一化,使得从x=1x=10积分结果为1。你是这样做的吗? - unutbu
1
我已经修改了我的答案,计算出离散概率质量函数。现在 pmf[0] = p(1) = 0.64525798 - unutbu
显示剩余3条评论

4

我不使用Python,所以我尝试描述解决方案的算法而不是冒险出现语法错误。这是一种暴力的离散反演方法。它应该很容易转换成Python代码。我假设数组采用基于0的索引。

设置:

  1. 生成一个大小为m的数组cdf,其中第一个条目为cdf [0] = 1,其余条目为cdf [i] = cdf [i-1] + 1 / (i +1) ** a

  2. 通过将每个条目除以cdf [m-1]来缩放所有条目--现在它们实际上是CDF值。

使用方法:

  • 通过生成均匀分布在[0,1]之间的随机数,并搜索cdf [],直到找到大于您的随机数的条目。返回索引+1作为您的x-值。

根据需要重复多次x-值。

例如,对于a,m = 2,10,我直接计算概率为:

[0.6452579827864142, 0.16131449569660355, 0.07169533142071269, 0.04032862392415089, 0.02581031931145657, 0.017923832855178172, 0.013168530260947229, 0.010082155981037722, 0.007966147935634743, 0.006452579827864143]

累积分布函数 (CDF) 是:

[0.6452579827864142, 0.8065724784830177, 0.8782678099037304, 0.9185964338278814, 0.944406753139338, 0.9623305859945162, 0.9754991162554634, 0.985581272236501, 0.9935474201721358, 1.0]

当生成时,如果我得到了一个均匀的结果为0.90,那么我将返回x=4,因为0.918...是大于我的均匀分布函数(CDF)的第一个条目。
如果您担心速度,可以构建一个别名表,但是由于几何衰减,线性搜索数组的早期终止概率非常高。例如,在给定的示例中,您将在第一个峰值上约有2/3的时间终止。

哎呀,我只用了两个小时(并阅读了你的答案)才意识到 OP 正在询问一个离散概率分布... - unutbu
这就是为什么我在询问是否可以采用范围区域来产生离散值。 - pjs

0
使用numpy.random.zipf函数,然后拒绝任何大于或等于m的样本。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接