从n维单位单形中随机均匀地采样

Question

从n维单位单形中随机均匀地采样是一种奇特的方式,可以说你需要n个随机数

• 他们都是非负面的,
• 他们总结为一个,和
• 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个减去它们的总和也会引入偏差.

维基百科提供了两种算法来正确地做到这一点:http: //en.wikipedia.org/wiki/Simplex#Random_sampling (虽然第二个目前声称只在实践中是正确的,不是理论上的.我希望能够清理它或者当我更好地理解这一点时澄清它.我最初在维基百科页面上发表了"警告:这样的论文声称以下是错误的",并且其他人将其变成了"仅在实践中作品"的警告.)

最后,问题是:您认为Mathematica中单纯形采样的最佳实现是什么(最好用经验证实它是正确的)？

相关问题

• 生成概率分布
• java随机百分比

Answer 1

这段代码可以工作:

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

基本上你只需要选择n - 1间隔[0,1]上的位置将其拆分,然后使用每个部分的大小Differences.

快速运行Timing表明它比Janus的第一个答案快一点.

Answer 2

经过一番挖掘,我发现这个页面提供了Dirichlet分布的很好的实现.从那里看来,遵循维基百科的方法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

Answer 3

我使用zdav:Dirichlet分布似乎是最简单的方法,并且zdav引用的Dirichlet分布的采样算法也出现在Dirichlet分布的维基百科页面上.

实现方式,首先进行完整的Dirichlet分配需要一些开销,因为你真正需要的只是n随机Gamma[1,1]样本.比较下面的
简单实现

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

完整Dirichlet实施

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}].

Answer 4

这是来自维基百科的第二个算法的简洁实现:

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 (最初它有联盟而不是排序@加入 - 后者稍微快一些.)

(有关证据证明这是正确的,请参阅评论!)