如果我想对截断整数幂律进行采样,我该使用Python中的哪个函数?
也就是说,给定两个参数a
和m
,在范围[1,m)
内生成一个随机整数x
,其遵循与1/x^a
成比例的分布。
我已经在numpy.random
中搜索过,但没有找到这个分布。
如果我想对截断整数幂律进行采样,我该使用Python中的哪个函数?
也就是说,给定两个参数a
和m
,在范围[1,m)
内生成一个随机整数x
,其遵循与1/x^a
成比例的分布。
我已经在numpy.random
中搜索过,但没有找到这个分布。
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()
1/x**a
。因此,没有在区间[1,2]、[2,3]等上进行积分。但是,我手动积分找到了const
和_ppf
的公式,即cdf
的反函数。我认为我做对了,但我可能错了。(我尝试了您的建议,但它将域移位到[1,11]
,所以如果我理解正确,这不符合基本的合理性检查。)顺便问一下,Spinal Tap在这里指的是什么? - unutbua,m = 2,10
进行了直接计算,p(1)应该为0.6452579827864142。这是您得到的结果吗? - pjs1/x**2
从x=1
到x=2
进行积分,并进行归一化,使得从x=1
到x=10
积分结果为1。你是这样做的吗? - unutbupmf[0] = p(1) = 0.64525798
。 - unutbu我不使用Python,所以我尝试描述解决方案的算法而不是冒险出现语法错误。这是一种暴力的离散反演方法。它应该很容易转换成Python代码。我假设数组采用基于0的索引。
设置:
生成一个大小为m
的数组cdf
,其中第一个条目为cdf [0] = 1
,其余条目为cdf [i] = cdf [i-1] + 1 / (i +1) ** a
。
通过将每个条目除以cdf [m-1]
来缩放所有条目--现在它们实际上是CDF值。
使用方法:
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]
x=4
,因为0.918...是大于我的均匀分布函数(CDF)的第一个条目。