稳定分布的随机数？

Question

如何在C＃中生成具有稳定分布的随机数？ Random类具有均匀分布。互联网上的许多其他代码显示正态分布。但我们需要稳定分布，即无限方差，也称为重尾分布。

原因是为了生成逼真的股票价格。在现实世界中，价格的巨大变化比在正态分布中更有可能发生。

有人知道将Random类输出转换为稳定分布的C＃代码吗？

编辑：嗯。确切的分布不如确保它将随机生成至少20 sigma这样巨大的sigma更重要。我们想测试交易策略对真正的重尾分布（这正是股市价格行为）的韧性。

我刚刚阅读了ZipFian和Cauchy的评论。由于必须选择一个，让我们选择Cauchy分布，但我也会尝试ZipFian进行比较。

- Wayne

---

我认为你需要找到所需精确分布的数学定义（正如@David所说，它可能应该有另一个名称），然后实现自己的函数来将均匀分布转换为所需分布。 - Zruty

在开始之前，您可以先学习一些基本的数学知识。Fisher变换（一种指标）进行高斯归一化 - 您要做的是相反的操作。基本上，您需要使用更大的随机空间（0-1000）来生成0-100的数字，并将更大的区域分配给距离0更远的数字。这实际上只是均匀分布的“高斯化”，也是我问那些愿意从事量化工作的程序员的问题之一 ;) 确保他们了解统计学101 ;) - TomTom

@DavidHeffernan：维基百科页面有些混乱。我认为正态分布是稳定分布的一个非常特殊的情况，因为它是唯一一个没有厚尾的稳定分布。然而，我已经很久很久没有学过统计学了，所以我可能在这里错了。 - Eric Lippert

@Wayne：你具体对哪种重尾分布感兴趣？比如说，柯西分布？ - Eric Lippert

@Eric 是的，我明白你的意思。事实上，只有在指定了特定的参数化时，这个问题才能真正得到回答。 - David Heffernan

Answer 1

通常的方法是：

• 选择一个稳定的、尾部较厚的分布，比如Cauchy分布。
• 查找所选分布的分位函数。

对于Cauchy分布而言，这将会是 p --> peak + scale * tan( pi * (p - 0.5) )。

现在你已经有了一种将均匀分布随机数转换为Cauchy分布随机数的方法。

明白了吗？具体细节请参考http://en.wikipedia.org/wiki/Inverse_transform_sampling。

注意：我学过统计学已经好久好久了。

更新：

我非常喜欢这个问题，刚刚写了一篇博客，详情请见：http://ericlippert.com/2012/02/21/generating-random-non-uniform-data/。

我的另一篇文章介绍了一些有趣的Zipfian分布示例，请参见：http://blogs.msdn.com/b/ericlippert/archive/2010/12/07/10100227.aspx。

- Eric Lippert

---

这在大多数分销版本中都很有效。至于正态分布，那就是另外一回事了！ - David Heffernan

Answer 2

如果您有兴趣使用Zipfian分布（通常用于建模来自科学或社会领域的过程），您可以按照以下步骤进行：

1. 选择分布的k值（偏斜度）
2. 预先计算累积分布的定义域（这只是一种优化）
3. 通过查找最接近的定义域值为分布生成随机值

示例代码：

List<int> domain = Enumerable.Range(0,1000);  // generate your domain
double skew  = 0.37; // select a skew appropriate to your domain
double sigma = domain.Aggregate(0.0d, (z,x) => x + 1.0 / Math.Pow(z+1, skew));
List<double> cummDist = domain.Select( 
      x => domain.Aggregate(0.0d, (z,y) => z + 1.0/Math.Pow(y, skew) * sigma));

现在，您可以通过从域中选择最接近的值来生成随机值：

Random rand = new Random();
double seek = rand.NextDouble();
int searchIndex = cummDist.BinarySearch(seek);
// return the index of the closest value from the distribution domain
return searchIndex < 0 ? (~searchIndex)-1 : searchIndex-1;

当然，您可以通过将实现分布域的逻辑从映射并返回该域值的过程中分离出来，从而概括整个过程。

- LBushkin

---

老兄，我很欣赏你优雅的代码，并且喜欢预先计算值的优化。你有没有想过如何将其转换为时间序列股票价格？ - Wayne

似乎柯西分布更适合于股票价格的随机数。你不同意吗？我会尝试两种方法，看哪一种效果更好。 - Wayne

@Wayne：我将把决定股价波动更准确的分布留给那些更有智慧和经验的人。关于金融市场建模有很多文章（维基百科只是冰山一角http://en.wikipedia.org/wiki/Modeling_and_analysis_of_financial_markets）。最终，根据您模拟中最重要的因素，只有您能够回答哪种分布最合适。 - LBushkin

Answer 3

下面的C#代码根据形状参数alpha和beta生成一个遵循稳定分布的随机数。我在Creative Commons Zero下将其发布到公共领域。

public static double StableDist(Random rand, double alpha, double beta){
    if(alpha<=0 || alpha>2 || beta<-1 || beta>1)
        throw new ArgumentException();
    var halfpi=Math.PI*0.5;
    var unif=NextDouble(rand);
    while(unif == 0.0)unif=NextDouble(rand);
    unif=(unif - 0.5) * Math.PI;
    // Cauchy special case
    if(alpha==1 && beta==0)
        return Math.Tan(unif);
    var expo=-Math.Log(1.0 - NextDouble(rand));
    var c=Math.Cos(unif);
    if(alpha == 1){
        var s=Math.Sin(unif);
        return 2.0*((unif*beta+halfpi)*s/c -
            beta * Math.Log(halfpi*expo*c/(
                unif*beta+halfpi)))/Math.PI;
    }
    var z=-Math.Tan(halfpi*alpha)*beta;
    var ug=unif+Math.Atan(-z)/alpha;
    var cpow=Math.Pow(c, -1.0 / alpha);
    return Math.Pow(1.0+z*z, 1.0 / (2*alpha))*
        (Math.Sin(alpha*ug)*cpow)*
        Math.Pow(Math.Cos(unif-alpha*ug)/expo, (1.0-alpha) / alpha);
}

private static double NextDouble(Random rand){
    // The default NextDouble implementation in .NET (see
    // https://github.com/dotnet/corert/blob/master/src/System.Private.CoreLib/shared/System/Random.cs)
    // is very problematic:
    // - It generates a random number 0 or greater and less than 2^31-1 in a 
    // way that very slightly biases 2^31-2.
    // - Then it divides that number by 2^31-1.
    // - The result is a number that uses roughly only 32 bits of pseudorandomness, 
    // even though `double` has 53 bits in its significand.
    // To alleviate some of these problems, this method generates a random 53-bit
    // random number and divides that by 2^53.  Although this doesn't fix the bias
    // mentioned above (for the default System.Random), this bias may be of
    // negligible importance for most purposes not involving security.
    long x=rand.Next(0,1<<30);
    x<<=23;
    x+=rand.Next(0,1<<23);
    return (double)x / (double)(1L<<53);
}

此外，我在另一篇文章中提供了稳定分布的伪代码。

- Peter O.

Answer 4

我手头有James Gentle的Springer专题书，随机数生成和蒙特卡罗方法，这是我统计学家妻子提供的。它在第105页讨论了稳定分布族：

“稳定分布族是一种灵活的、通常具有重尾分布的分布族。这个族群包括正态分布在一个参数的极端值和柯西分布在另一个极端值。 Chambers、Mallows和Stuck（1976）提供了一种从稳定分布中生成偏差的方法。（注意辅助函数D2中常数的错误，用于评估（e^x-1）/x）。他们的方法在IMSL库中使用。对于对称稳定分布，Devroye（1986）指出可以通过利用对称稳定与Fejer-de la Vallee Poissin分布的关系来开发更快的方法。 Buckle（1995）展示了如何在数据条件下模拟稳定分布的参数。”

从通用稳定分布中生成偏差很困难。如果您需要这样做，我建议使用IMSL等库。我不建议您自己尝试。

然而，如果您正在寻找稳定族中的特定分布，例如柯西分布，则可以使用Eric描述的方法，即概率积分变换。只要您能够以闭合形式写出分布函数的反函数，就可以使用此方法。

- David Heffernan