从 n 维单位单纯形中随机均匀采样答案

【问题标题】：Sample uniformly at random from an n-dimensional unit simplex从 n 维单位单纯形中随机均匀采样
【发布时间】：2011-03-01 22:47:57
【问题描述】：

从一个 n 维单位单纯形中均匀随机抽样是说你想要 n 个随机数这样的奇特方式

它们都是非负的，
它们总和为 1，并且
n 个非负数总和为 1 的每个可能向量的可能性均等。

在 n=2 的情况下，您希望从位于正象限的线 x+y=1（即 y=1-x）的线段中均匀采样。在 n=3 的情况下，您从平面 x+y+z=1 的三角形部分进行采样，该部分位于 R3 的正八分圆中：

（图片来自http://en.wikipedia.org/wiki/Simplex。）

请注意，选择 n 个统一的随机数，然后将它们归一化以使它们总和为 1 是行不通的。你最终会偏向不太极端的数字。

类似地，选取 n-1 个均匀随机数，然后将第 n 个减去它们的总和后取第 n 个，也会引入偏差。

Wikipedia 提供了两种算法来正确执行此操作：http://en.wikipedia.org/wiki/Simplex#Random_sampling （虽然第二个目前声称仅在实践中是正确的，而不是在理论上。我希望在我更好地理解这一点时对其进行清理或澄清。我最初陷入了“警告：某某纸声称该维基百科页面上的以下内容是错误的”，其他人将其变成了“仅在实践中有效”的警告。）

最后，问题：您认为 Mathematica 中单纯形采样的最佳实现方式是什么（最好通过经验确认它是正确的）？

相关问题

【问题讨论】：

似乎有几种方法可以正常工作 - 唯一真正的区别在于速度和可读性。除了“最佳”之外，您的标准是什么？
速度和可读性是很好的标准！简洁可能是另一个。如果您的实现有任何用途，请继续并将其作为答案发布。
我认为维基百科的警告有点虚假；引用的论文的作者担心这个问题的离散化版本的完美一致性。从数学的角度来看，所描述的第二种算法是完全正确的，如果您准备将“来自 [0, 1] 的随机浮点数”视为“随机实数”的足够好的近似值，那么在实践中应该可以很好地工作[0, 1]'中的数字。
采样链接失效
Related.

标签： math random wolfram-mathematica

【解决方案1】：

老问题，我迟到了，但如果有效实施，这种方法比公认的答案要快得多。

在 Mathematica 代码中： #/总计[#,{2}]&@Log@RandomReal[{0,1},{n,d}]

用简单的英语，您生成 n 行 * d 列随机分布在 0 和 1 之间。然后取所有内容的 Log。然后对每一行进行归一化，将行中的每个元素除以行总数。现在你有 n 个样本均匀分布在 (d-1) 维单纯形上。

如果在这里找到这个方法：https://mathematica.stackexchange.com/questions/33652/uniformly-distributed-n-dimensional-probability-vectors-over-a-simplex

我承认，我不确定它为什么有效，但它通过了我能想到的所有统计测试。如果有人能证明为什么这种方法有效，我很乐意看到！

【讨论】：

【解决方案2】：

我创建了一个算法，用于在单纯形上进行均匀随机生成。您可以在以下链接中找到论文中的详细信息： http://www.tandfonline.com/doi/abs/10.1080/03610918.2010.551012#.U5q7inJdVNY

简单地说，您可以使用以下递归公式来找到 n 维单纯形上的随机点：

x₁=1-R₁^1/n-1

x_k=(1-Σ_i=1^{k i>}x_i)(1-R_{k i>}^1/nk), k=2, ..., n-1

x_n=1-Σ_i=1^n-1x_i

其中 R_i 是 0 到 1 之间的随机数。

现在我正在尝试制作一种算法来从受约束的单纯形中生成随机均匀样本。这是单纯形和凸体之间的交集。

【讨论】：

你知道你会如何改变这个算法来增加一个维度的权重吗？

【解决方案3】：

这段代码可以工作：

samples[n_] := Differences[Join[{0}, Sort[RandomReal[Range[0, 1], n - 1]], {1}]]

基本上，您只需在间隔[0,1] 上选择n - 1 位置将其拆分，然后使用Differences 获取每个部分的大小。

Timing 对此的快速运行表明它比 Janus 的第一个答案要快一些。

【讨论】：

谢谢！我认为这与我发布的几乎同构。感谢您的速度比较，顺便说一句！
哦，确实如此！不知怎的，我以为你的不一样，但现在我再看了一遍，看起来还是一样。
是否有关于非单位单纯形甚至（单纯）凸多面体的概括？
@Orient 你不能把结果放大到你想要的任何大小吗？
@BenAlpert 我可以简单地将结果放大到，一方面，根据产生的单纯形的超体积对所需凸多面体进行三角剖分，但另一方面，我不能简单地变形仅通过线性变换将单位单纯形转化为任意单纯形，并通过这样做来保持空间分布的均匀性。

【解决方案4】：

这是来自Wikipedia 的第二个算法的简洁实现：

SimplexSample[n_] := Rest@# - Most@# &[Sort@Join[{0,1}, RandomReal[{0,1}, n-1]]]

改编自这里：http://www.mofeel.net/1164-comp-soft-sys-math-mathematica/14968.aspx （最初它使用 Union 而不是 Sort@Join ——后者稍微快一些。）

（请参阅 cmets 以获得一些证明这是正确的证据！）

【讨论】：

我刚刚进行了测试，它似乎工作正常 - 值非常统一并且都在正确的行上。我不知道为什么它被否决了。
投反对票是我的，是偶然的；我道歉。如果您可以对您的答案进行简单的编辑，我会投赞成票。这种方法在我看来是正确的。
我只是要发布同样的方法。在我的机器上，它比 Gamma[1,1] 采样方法快八倍。
不用Rest - Most，也可以直接调用内置的Differences，也是一样的。

【解决方案5】：

经过一番挖掘，我找到了this page，它很好地实现了狄利克雷分布。从那里开始，遵循维基百科的方法 1 似乎非常简单。这似乎是最好的方法。

作为测试：

In[14]:= RandomReal[DirichletDistribution[{1,1}],WorkingPrecision->25]
Out[14]= {0.8428995243540368880268079,0.1571004756459631119731921}
In[15]:= Total[%]
Out[15]= 1.000000000000000000000000

100 个样本的图：

alt text http://www.public.iastate.edu/~zdavkeos/simplex-sample.png

【讨论】：

【解决方案6】：

我支持 zdav：Dirichlet 分布似乎是最简单的方法，zdav 所指的对 Dirichlet 分布进行采样的算法也在 Dirichlet distribution 的维基百科页面上提供。

在实现方面，首先进行完整的狄利克雷分布有点开销，因为您真正需要的只是n随机Gamma[1,1]样本。下面比较
简单实现

SimplexSample[n_, opts:OptionsPattern[RandomReal]] :=
  (#/Total[#])& @ RandomReal[GammaDistribution[1,1],n,opts]

完整的狄利克雷实现

DirichletDistribution/:Random`DistributionVector[
 DirichletDistribution[alpha_?(VectorQ[#,Positive]&)],n_Integer,prec_?Positive]:=
    Block[{gammas}, gammas = 
        Map[RandomReal[GammaDistribution[#,1],n,WorkingPrecision->prec]&,alpha];
      Transpose[gammas]/Total[gammas]]

SimplexSample2[n_, opts:OptionsPattern[RandomReal]] := 
  (#/Total[#])& @ RandomReal[DirichletDistribution[ConstantArray[1,{n}]],opts]

时机

Timing[Table[SimplexSample[10,WorkingPrecision-> 20],{10000}];]
Timing[Table[SimplexSample2[10,WorkingPrecision-> 20],{10000}];]
Out[159]= {1.30249,Null}
Out[160]= {3.52216,Null}

所以完整的 Dirichlet 慢了 3 倍。如果您一次需要 m>1 个采样点，您可能会通过 (#/Total[#]&)/@RandomReal[GammaDistribution[1,1],{m,n}] 获得更多胜利。

【讨论】：