产生幂律分布的随机数发生器?
我正在为C++命令行Linux应用程序编写一些测试。我想用幂律/长尾分布生成一堆整数。意思是,我经常收到一些数字,但其中大多数并不常见。产生幂律分布的随机数发生器?
理想情况下,我只能用rand()或stdlib随机函数之一来使用一些魔术方程。如果没有,一个简单易用的C/C++将会非常棒。
谢谢!
这个page at Wolfram MathWorld讨论了如何从均匀分布(这是大多数随机数发生器提供的)得到幂律分布。
简短的回答(在上面的链接派生):
x = [(x1^(n+1) - x0^(n+1))*y + x0^(n+1)]^(1/(n+1))
其中ý是均匀的变量,Ñ是分布功率,X0和X1限定的范围内分布,而x是你的幂律分布变量。
如果您知道您想要的分布(称为概率分布函数(PDF))并将其正确化,则可以将其整合以获得累积分布函数(CDF),然后将CDF(如果可能)从统一的[0,1]
分配到你想要的转换。
所以,你首先定义你想要的发行版。
P = F(x)
(对于x在[0,1]),然后积分得到
C(y) = \int_0^y F(x) dx
如果能倒你
y = F^{-1}(C)
于是呼rand()
和堵塞结果作为C
在最后一行并使用y。
这个结果被称为抽样的基本定理。这是一个麻烦,因为规范化要求和需要分析反转功能。
或者,您可以使用拒绝技术:在所需范围内均匀地抛出一个数字,然后抛出另一个数字并与第一次抛出的位置处的PDF进行比较。如果第二次投掷超过PDF,则拒绝。对于具有很多低概率区域的PDF,倾向于效率低下,如长尾巴的那些......
中间方法涉及通过强力反转CDF:将CDF作为查找表存储,并执行反转查找以获得结果。
这里真正臭气熏天就是这么简单x^-n
分布都在范围[0,1]
非normalizable,所以你不能使用采样定理。尝试(x + 1)^ - n改为...
我无法评论产生幂律分布所需的数学(其他职位有建议),但我建议您熟悉<random>
中的TR1 C++标准库随机数设施。这些提供比std::rand
和std::srand
更多的功能。新系统为发电机,发动机和配电系统指定了一个模块化API,并提供了一堆预设。
所包含的分配预设有:
uniform_int
bernoulli_distribution
geometric_distribution
poisson_distribution
binomial_distribution
uniform_real
exponential_distribution
normal_distribution
gamma_distribution
当你定义你的幂律分布,你应该能够与现有的发电机和引擎插入。本书Pete Becker的C++标准库扩展对<random>
有很大的帮助。
Here is an article有关如何创建其他分布(与柯西,卡方,学生吨和费雪˚F例子)
我只是想进行实际模拟作为补充,(理所当然)接受的答案。尽管在R中,代码非常简单,可以成为(伪) - 伪代码。在接受的答案和其他的Wolfram MathWorld formula之间
一个微小的差异,也许更常见的,方程是一个事实,即幂指数n
(通常称为阿尔法)不进行明确的负号。所以选择的alpha值必须是负数,通常在2和3之间。
x0
和x1
表示分布的上限和下限。
所以在这里,它是:
x1 = 5 # Maximum value
x0 = 0.1 # It can't be zero; otherwise X^0^(neg) is 1/0.
alpha = -2.5 # It has to be negative.
y = runif(1e5) # Number of samples
x = ((x1^(alpha+1) - x0^(alpha+1))*y + x0^(alpha+1))^(1/(alpha+1))
hist(x, prob = T, breaks=40, ylim=c(0,10), xlim=c(0,1.2), border=F,
col="yellowgreen", main="Power law density")
lines(density(x), col="chocolate", lwd=1)
lines(density(x, adjust=2), lty="dotted", col="darkblue", lwd=2)
或对数刻度绘制:
h = hist(x, prob=T, breaks=40, plot=F)
plot(h$count, log="xy", type='l', lwd=1, lend=2,
xlab="", ylab="", main="Density in logarithmic scale")
下面是数据汇总:
> summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.1000 0.1208 0.1584 0.2590 0.2511 4.9388
当极限值为0和无限时,这是否工作? – Peaceful 2015-01-06 06:11:00
小额外细节:** y **是[0,1]范围内的均匀变量。 – 2017-01-12 03:22:01