产生幂律分布的随机数发生器？

Question

我正在为C++命令行Linux应用程序编写一些测试。我想用幂律/长尾分布生成一堆整数。意思是，我经常收到一些数字，但其中大多数并不常见。

理想情况下，我只能用rand（）或stdlib随机函数之一来使用一些魔术方程。如果没有，一个简单易用的C/C++将会非常棒。

谢谢！

Answer 1

这个page at Wolfram MathWorld讨论了如何从均匀分布（这是大多数随机数发生器提供的）得到幂律分布。

简短的回答（在上面的链接派生）：

x = [(x1^(n+1) - x0^(n+1))*y + x0^(n+1)]^(1/(n+1)) 

其中ý是均匀的变量，Ñ是分布功率，X0和X1限定的范围内分布，而x是你的幂律分布变量。

Answer 2

如果您知道您想要的分布（称为概率分布函数（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改为...

Answer 3

我无法评论产生幂律分布所需的数学（其他职位有建议），但我建议您熟悉<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例子）

Answer 4

我只是想进行实际模拟作为补充，（理所当然）接受的答案。尽管在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